import numpy as np
from flopy import version as vs
if vs.__version__ >= '3.2.4':
    from flopy.pakbase import Package
if vs.__version__ == '3.2.3':
    from flopy.mbase import Package


class ModflowFmp(Package):
    """
    'Farm Process Package Class'

    Parameters
    ----------
    model : model object
        The model object (of type :class:'flopy.modflow.mf.Modflow') to which
        this package will be added.
    mnwmax : integer
        Maximum number of multi-node wells (MNW) to be simulated (default is 0)
    ipakcb : integer
        A flag that is used to determine if cell-by-cell budget data should be saved
            if ipakcb > 0, then it is the unit number to which MNW cell-by-cell flow terms will be
                recorded whenever cell-by-cell budget data are written to a file (as determined by the
                output control options of MODFLOW).
            if ipakcb = 0, then MNW cell-by-cell flow terms will not be printed or recorded.
            if ipakcb < 0, then well injection or withdrawal rates and water levels in the well and
                its multiple cells will be printed in the main MODFLOW listing (output) file whenever cell-by-cell
                budget data are written to a file (as determined by the output control options of MODFLOW).
        (default is 0)
    mnwprnt : integer
        Flag controlling the level of detail of information about multi-node wells to be written to the
        main MODFLOW listing (output) file. If MNWPRNT = 0, then only basic well information will be
        printed in the main MODFLOW output file; increasing the value of MNWPRNT yields more information,
        up to a maximum level of detail corresponding with MNWPRNT = 2.
        (default is 0)
    aux : list of strings
    (note: not sure if the words AUX or AUXILIARY are necessary)
        in the style of AUXILIARY abc or AUX abc where
        abc is the name of an auxiliary parameter to be read for each multi-node well as part
        of dataset 4a. Up to five parameters can be specified, each of which must be preceded
        by AUXILIARY or AUX. These parameters will not be used by the MNW2 Package, but
        they will be available for use by other packages.
        (default is None)
    wellid : array of strings (shape = (MNWMAX))
        The name of the wells. This is a unique identification label for each well.
        The text string is limited to 20 alphanumeric characters. If the name of the well includes
        spaces, then enclose the name in quotes.
        (default is None)
    nnodes : integers
        The number of cells (nodes) associated with this well. NNODES normally is > 0, but for the
        case of a vertical borehole, setting NNODES < 0 will allow the user to specify the elevations
        of the tops and bottoms of well screens or open intervals (rather than grid layer numbers),
        and the absolute value of NNODES equals the number of open intervals (or well screens) to be
        specified in dataset 2d. If this option is used, then the model will compute the layers in which
        the open intervals occur, the lengths of the open intervals, and the relative vertical position
        of the open interval within a model layer (for example, see figure 14 and related discussion)
        (default is None)
    losstype : string
        The user-specified model for well loss. The following loss types are currently supported.
            NONE there are no well corrections and the head in the well is assumed to equal the head
                in the cell. This option (hWELL = hn) is only valid for a single-node well (NNODES = 1).
                (This is equivalent to using the original WEL Package of MODFLOW, but specifying the single-node
                well within the MNW2 Package enables the use of constraints.)
            THIEM this option allows for only the cell-to-well correction at the well based on the Thiem (1906)
                equation; head in the well is determined from equation 2 as (hWELL = hn + AQn), and the model
                computes A on the basis of the user-specified well radius (Rw) and previously defined values of
                cell transmissivity and grid spacing. Coefficients B and C in equation 2 are automatically
                set = 0.0. User must define Rw in dataset 2c or 2d.
            SKIN this option allows for formation damage or skin corrections at the well:
                hWELL = hn + AQn + BQn (from equation 2), where A is determined by the model from the value of
                Rw, and B is determined by the model from Rskin and Kskin. User must define Rw, Rskin, and Kskin in
                dataset 2c or 2d.
            GENERAL head loss is defined with coefficients A, B, and C and power exponent P
                (hWELL = hn + AQn + BQn + CQnP). A is determined by the model from the value of Rw.
                must define Rw, B, C, and P in dataset 2c or 2d. A value of P = 2.0 is suggested if no other
                data are available (the model allows 1.0 <= P <= 3.5). Entering a value of C = 0 will result
                in a linear model in which the value of B is entered directly (rather than entering properties
                of the skin, as with the SKIN option).
            SPECIFYcwc the user specifies an effective conductance value (equivalent to the combined
                effects of the A, B, and C well-loss coefficients expressed in equation 15) between the well and
                the cell representing the aquifer, CWC. User must define CWC in dataset 2c or 2d. If there are
                multiple screens within the grid cell or if partial penetration corrections are to be made, then
                the effective value of CWC for the node may be further adjusted automatically by MNW2.
                (default is None)
    pumploc : integer
        The location along the borehole of the pump intake (if any).
        If PUMPLOC = 0, then either there is no pump or the intake location (or discharge point for an
        injection well) is assumed to occur above the first active node associated with the multi-node well
        (that is, the node closest to the land surface or to the wellhead). If PUMPLOC > 0, then the cell in
        which the intake (or outflow) is located will be specified in dataset 2e as a LAY-ROW-COL grid location.
        For a vertical well only, specifying PUMPLOC < 0, will enable the option to define the vertical position of
        the pump intake (or outflow) as an elevation in dataset 2e (for the given spatial grid location [ROW-COL]
        defined for this well in 2d).
        (default is 0)
    qlimit : integer
        Indicates whether the water level (head) in the well will be used to constrain
        the pumping rate. If Qlimit = 0, then there are no constraints for this well. If Qlimit > 0, then
        pumpage will be limited (constrained) by the water level in the well, and relevant parameters are
        constant in time and defined below in dataset 2f. If Qlimit < 0, then pumpage will be limited
        (constrained) by the water level in the well, and relevant parameters can vary with time and are
        defined for every stress period in dataset 4b.
        (default is 0)
    ppflag : integer
        Flag that determines whether the calculated head in the well will be corrected for the
        effect of partial penetration of the well screen in the cell. If PPFLAG = 0, then the head in the
        well will not be adjusted for the effects of partial penetration. If PPFLAG > 0, then the head in
        the well will be adjusted for the effects of partial penetration if the section of well containing
        the well screen is vertical (as indicated by identical row-column locations in the grid). If
        NNODES < 0 (that is, the open intervals of the well are defined by top and bottom elevations),
        then the model will automatically calculate the fraction of penetration for each node and the
        relative vertical position of the well screen. If NNODES > 0, then the fraction of penetration for
        each node must be defined in dataset 2d (see below) and the well screen will be assumed to be
        centered vertically within the thickness of the cell (except if the well is located in the uppermost
        model layer that is under unconfined conditions, in which case the bottom of the well screen will be
        assumed to be aligned with the bottom boundary of the cell and the assumed length of well screen will
        be based on the initial head in that cell).
        (default is 0)
    pumpcap : integer
        A flag and value that determines whether the discharge of a pumping (withdrawal) well (Q < 0.0)
        will be adjusted for changes in the lift (or total dynamic head) with time. If PUMPCAP = 0,
        then the discharge from the well will not be adjusted on the basis of changes in lift. If PUMPCAP > 0
        for a withdrawal well, then the discharge from the well will be adjusted on the basis of the lift, as
        calculated from the most recent water level in the well. In this case, data describing the head-capacity
        relation for the pump must be listed in datasets 2g and 2h, and the use of that relation can be switched
        on or off for each stress period using a flag in dataset 4a. The number of entries (lines) in dataset 2h
        corresponds to the value of PUMPCAP. If PUMPCAP does not equal 0, it must be set to an integer value of
        between 1 and 25, inclusive.
        (default is 0)
    lay_row_col : list of arrays (shape = (NNODES,3), length = MNWMAX)
        Layer, row, and column numbers of each model cell (node) for the current well. If NNODES > 0,
        then a total of NNODES model cells (nodes) must be specified for each well (and dataset 2d must
        contain NNODES records). In the list of nodes defining the multi-node well, the data list must be
        constructed and ordered so that the first node listed represents the node closest to the wellhead,
        the last node listed represents the node furthest from the wellhead, and all nodes are listed in
        sequential order from the top to the bottom of the well (corresponding to the order of first to
        last well nodes). A particular node in the grid can be associated with more than one multi-node well.
        (default is None)
    ztop_zbotm_row_col : list of arrays (shape = (abs(NNODES),2), length = MNWMAX)
        The top and bottom elevations of the open intervals (or screened intervals) of a vertical well.
        These values are only read if NNODES < 0 in dataset 2a. The absolute value of NNODES indicates
        how many open intervals are to be defined, and so must correspond exactly to the number of records
        in dataset 2d for this well. In the list of intervals defining the multi-node well, the data list
        must be constructed and ordered so that the first interval listed represents the shallowest one,
        the last interval listed represents the deepest one, and all intervals are listed in sequential
        order from the top to the bottom of the well. If an interval partially or fully intersects a model
        layer, then a node will be defined in that cell. If more than one open interval intersects a
        particular layer, then a length-weighted average of the cell-to-well conductances will be used to
        define the well-node characteristics; for purposes of calculating effects of partial penetration,
        the cumulative length of well screens will be assumed to be centered vertically within the thickness
        of the cell. If the well is a single-node well by definition of LOSSTYPE = NONE and the defined open
        interval straddles more than one model layer, then the well will be associated with the cell where
        the center of the open interval exists.
        (default is None)

        if losstype != None (see losstype for definitions)
            rw : float
                (default is 0)
            rskin : float
                (default is 0)
            kskin : float
                (default is 0)
            b : float
                (default is 0)
            c : float
                (default is 0)
            p : float
                (default is 0)
            cwc  float
                (default is 0)
    pp : float
        the fraction of partial penetration for this cell
        (see PPFLAG in dataset 2b). Only specify if PPFLAG > 0 and NNODES > 0.
        (default is 1)
    itmp : integer
        For reusing or reading multi-node well data; it can change each stress period.
        ITMP must be >= 0 for the first stress period of a simulation.
        if ITMP > 0, then ITMP is the total number of active multi-node wells simulated during the stress period,
        and only wells listed in dataset 4a will be active during the stress period. Characteristics of each
        well are defined in datasets 2 and 4.
        if ITMP = 0, then no multi-node wells are active for the stress period and the following dataset is
        skipped.
        if ITMP < 0, then the same number of wells and well information will be reused from the previous stress
        period and dataset 4 is skipped.
        (default is 0)
    wellid_qdes : list of arrays (shape = (NPER,MNWMAX,2))
        the actual (or maximum desired, if constraints are to be applied) volumetric pumping rate
        (negative for withdrawal or positive for injection) at the well (L3/T). Qdes should be set to 0
        for nonpumping wells. If constraints are applied, then the calculated volumetric withdrawal or
        injection rate may be adjusted to range from 0 to Qdes and is not allowed to switch directions
        between withdrawal and injection conditions during any stress period. When PUMPCAP > 0, in the
        first stress period in which Qdes is specified with a negative value, Qdes represents the maximum
        operating discharge for the pump; in subsequent stress periods, any different negative values of
        Qdes are ignored, although values are subject to adjustment for CapMult. If Qdes >= 0.0, then
        pump-capacity adjustments are not applied.
        (default is None)
    extension : string
        Filename extension (default is 'fmp')
    unitnumber : int
        File unit number (default is 1001).

    Attributes
    ----------

    Methods
    -------

    See Also
    --------

    Notes
    -----


    """

    def __init__(self, model, npfwl = 0, mxl = 0, mxactw=0, nfarms=0, ncrops=0, nsoils=0,
                 ifrmfl = 1, irtfl = 3, icufl = -1, ipfl = 2, iftefl = 1, iieswfl = 1, ieffl = 1,
                 iebfl = 0, irotfl = 0, ideffl = -2, iben = 1, icost = 1, iallotgw = 1,
                 iccfl = 1, inrdfl = 0, mxnrdt = 0, isrdfl = 1, irdfl = 0, isrrfl = 0, irrfl = -1, iallotsw = 0, pclose = 0.1, ifwlcb = 1, ifnrcb = 1, isdpfl = 1, ifbpfl = 1, ietpfl = 1, irtpfl = 2, iopfl = 2, ipapfl = 1, itmp = None,
                 farmwells = None, gse = None, fid = None, ofe = None, farmcost = None, nonrouted = None, allotsw = None, allotgw = None, srswlocation = None,
                     sid = None, soillist = None, cid = None, root = None, cropevap = None, cropfies = None,
                 psi = None, cropgrowth = None, cropirr = None, kc = None, etref = None, precip = None, ifallow = None, cropcost = None, climate = None,
                 extension='fmp', unitnumber=1001):
        """
        Package constructor
        """
        Package.__init__(self, model, extension, 'FMP',
                         unitnumber)  # Call ancestor's init to set self.parent, extension, name, and unit number

        self.url = 'fmp.htm'
        self.model = model
        self.nrow, self.ncol, self.nlay, self.nper = self.parent.nrow_ncol_nlay_nper

        self.heading = '# Farm process (FMP) file for MODFLOW, generated by Flopy'
        self.npfwl = npfwl
        self.mxl = mxl
        self.mxactw = mxactw
        self.nfarms = nfarms
        self.ncrops = ncrops
        self.nsoils = nsoils

        self.ifrmfl = ifrmfl
        self.irtfl = irtfl
        self.icufl = icufl
        self.ipfl = ipfl
        self.iftefl = iftefl
        self.iieswfl = iieswfl
        self.ieffl = ieffl

        self.iebfl = iebfl
        self.irotfl = irotfl
        self.ideffl = ideffl
        self.iben = iben
        self.icost = icost
        self.iallotgw = iallotgw

        self.iccfl = iccfl

        self.inrdfl = inrdfl
        self.mxnrdt = mxnrdt
        self.isrdfl = isrdfl
        self.irdfl = irdfl
        self.isrrfl = isrrfl
        self.irrfl = irrfl
        self.iallotsw = iallotsw
        self.pclose = pclose

        self.ifwlcb = ifwlcb
        self.ifnrcb = ifnrcb
        self.isdpfl = isdpfl
        self.ifbpfl = ifbpfl
        self.ietpfl = ietpfl
        self.irtpfl = irtpfl
        self.iopfl = iopfl
        self.ipapfl = ipapfl

        self.itmp = itmp

        self.farmwells = farmwells
        self.gse = gse
        self.fid = fid

        self.ofe = ofe
        self.farmcost = farmcost
        self.nonrouted = nonrouted
        self.allotgw = allotgw
        self.allotsw = allotsw
        self.srswlocation = srswlocation

        self.sid = sid
        self.soillist = soillist
        self.cid = cid
        self.root = root
        self.cropevap = cropevap
        self.cropfies = cropfies
        self.psi = psi
        self.cropgrowth = cropgrowth
        self.cropirr = cropirr
        self.kc = kc
        self.etref = etref
        self.precip = precip
        self.ifallow = ifallow
        self.cropcost = cropcost

        self.climate = climate

        # Ci vuole un metodo che scriva i file esterni .in
        # nella sotto cartella data_fmp/

        self.parent.add_package(self)

    def write_in_files(self):
        """
        Write the external files (*.in) that are called by file input (*.fmp) writer.

        Returns
        -------
        None

        """
        import os
        datapath = self.model.model_ws + '\\DATA_FMP'
        if not os.access(datapath, os.F_OK):
            os.mkdir(datapath)

        datapath = datapath + '\\'

        # GSE - ground surface
        fgse = open(datapath + 'GSE.in', 'w')
        for nr in range(self.nrow):
            for nc in range(self.ncol):
                fgse.write('{0:4e}  '.format(float(self.gse[nr,nc])))
            fgse.write('\n')
        fgse.close()
        # #
        # # FID - farms id
        ffid = open(datapath + 'FID.in', 'w')
        for nr in range(self.nrow):
            for nc in range(self.ncol):
                ffid.write('{0:4d}'.format(int(self.fid[nr,nc])))   #GDF In the original file, the format was {0:2d} (see issue #197)
            ffid.write('\n')
        ffid.close()
        # #
        # OFE   - for each farm OFE - Item 7
        fofe = open(datapath + 'OFE.in', 'w')
        for i in self.ofe.keys() :
            fofe.write('%d  '%i)
            for j in self.root.keys():
                fofe.write('%10f '% self.ofe[i]  )
            fofe.write('\n')
        fofe.close()
        # #
        # FARM COST - for each farm [Farm-ID GWcost1 GWcost2 GWcost3 GWcost4 SWcost1 SWcost2 SWcost3 SWcost4] - Item 21
        ffcost = open(datapath + 'FARMCOST.in', 'w')
        for i in self.farmcost.keys():
            ffcost.write('%4d  %10f  %10f  %10f  %10f  %10f  %10f  %10f  %10f \n'%(i, self.farmcost[i][0], self.farmcost[i][1], self.farmcost[i][2], self.farmcost[i][3], self.farmcost[i][4], self.farmcost[i][5], self.farmcost[i][6], self.farmcost[i][7] ))
        ffcost.close()
        # #
        # #
        # Semi-Routed Locations - for each farm [Farm-ID Row Column Segment Reach]  - Item 21
        if self.srswlocation is not None:
            fsrloc = open(datapath + 'SRSWLOC.in', 'w')
            for i in self.srswlocation.keys():
                fsrloc.write('%4d  %4d  %4d  %4d  %4d \n'%(i, self.srswlocation[i][0], self.srswlocation[i][1], self.srswlocation[i][2], self.srswlocation[i][3]  ))
            fsrloc.close()
        #


        # # SID - soils id
        fsid = open(datapath + 'SID.in', 'w')
        for nr in range(self.nrow):
            for nc in range(self.ncol):
                fsid.write('{0:4d}'.format(int(self.sid[nr,nc])))   #GDF In the original file, the format was {0:2d} (see issue #197)
            fsid.write('\n')
        fsid.close()
        # #
        # SOILLIST - soils list
        fslist = open(datapath + 'SOILLIST.in', 'w')
        for i in self.soillist.keys():
            fslist.write('{0:3d}{1:10f}  {2:15s} \n'.format(self.soillist[i][0], self.soillist[i][1],self.soillist[i][2]))
        fslist.close()
        #
        # # CID - crops id
        fcid = open(datapath + 'CID.in', 'w')
        for nr in range(self.nrow):
            for nc in range(self.ncol):
                fcid.write('{0:4d}'.format(int(self.cid[nr,nc])))   #GDF In the original file, the format was {0:2d} (see issue #197)
            fcid.write('\n')
        fcid.close()
        # #
        # CROPEVAP   - for each crop [FTR, FEP, FEI]
        fevap = open(datapath + 'CROPEVAP.in', 'w')
        for i in self.cropevap.keys():
            fevap.write('%4d  %10f  %10f  %10f \n'%(i, self.cropevap[i][0], self.cropevap[i][1], self.cropevap[i][2]  ))
        fevap.close()
        # #
        # CROP Fraction of in-efficient losses   -   [FIESWP, FIESWI]
        ffies = open(datapath + 'CROPFIES.in', 'w')
        for i in self.cropfies.keys():
            ffies.write('%4d  %4f  %4f \n'%(i, self.cropfies[i][0], self.cropfies[i][1]  ))
        ffies.close()
        # #
        # Consumptive Use Concept variables CROP -   [PSI1 to PSI4] - Item 14
        fpsi = open(datapath + 'PSI.in', 'w')
        for i in self.psi.keys():
            fpsi.write('%4d  %4f  %4f  %4f  %4f  \n'%(i, self.psi[i][0], self.psi[i][1], self.psi[i][2], self.psi[i][3]  ))   #GDF In the original file, the format for self.psi[i][0] was %4d (see issue #197) 
        fpsi.close()
        # #
        if self.irtfl == 2:
            # Root depth specified for each stress period.
            # # For each stress period: CropId Root  - Item 29
            #froot = open(datapath + 'ROOT.in', 'w')   #GDF (see issue #194)
            for sp in self.kc.keys():
                froot = open(datapath + 'ROOT_%s.in'%str(sp), 'w')   #GDF (see issue #194)
                for i, cp in enumerate(self.cropirr.keys()):
                    froot.write('%4d %4f\n'%(cp, self.root[sp][i] ))
                froot.close()   #GDF (see issue #194)
        # #
        if self.icufl == -1:
            # Consumptive Use Concept variables CROP
            # Kc read for each crop for stress period
            # ET array for each stress period
            #fkc = open(datapath + 'KC.in', 'w')   #GDF (see issue #194)
            #print 'sp in kc ', self.kc.keys()     #GDF (see issue #194)
            for sp in self.kc.keys():
                # Reference ET as array (nrow,ncol) for each SP, saved in a separate file for each SP
                fetref = open(datapath + 'ET_%s.in'%str(sp), 'w')
                for nr in range(self.nrow):
                    for nc in range(self.ncol):
                        #fetref.write('{0:4e}  '.format(float(self.etref[sp][nr,nc])))
                        fetref.write('%4.2e  '%self.etref[sp][nr,nc])
                    fetref.write('\n')
                fetref.close()

                # For each stress period: CropId Kc NONIRR  - Item 30a
                fkc = open(datapath + 'KC_%s.in'%str(sp), 'w')   #GDF (see issue #194)
                for i, cp in enumerate(self.cropirr.keys()):
                    fkc.write('%4d %4f %4d \n'%(cp, self.kc[sp][i], self.cropirr[cp][0]))

                fkc.close()   #GDF (see issue #194)

        if self.icufl == 3:
            # Consumptive Use Concept variables CROP -   [BaseT, MinCutT ,MaxCutT ,C0 ,C1 ,C2, C3, BegRootD, MaxRootD, RootGC, NONIRR,] - Item 15
            fgrowth = open(datapath + 'GROWTH.in', 'w')
            for i in self.cropgrowth.keys():
                fgrowth.write('%4d  %4f  %4f  %4f  %4f  %4f  %4f  %4f  %4f  %4f  %4f  %4d  \n'%(i, self.cropgrowth[i][0], self.cropgrowth[i][1], self.cropgrowth[i][2], self.cropgrowth[i][3], self.cropgrowth[i][4], self.cropgrowth[i][5], self.cropgrowth[i][6], self.cropgrowth[i][7], self.cropgrowth[i][8], self.cropgrowth[i][9], self.cropgrowth[i][10]  ))
            fgrowth.close()
        # #
        # Climate Data  - [TimeSeriesStep MaxT MinT Precip ETref]  - Item 16
        if self.icufl == 3:
            fclimate = open(datapath + 'CLIMATE.in', 'w')
            for i in range(len(self.climate)):
                fclimate.write('%4.2e  %4.2e  %4.2e  %4.2e  %4.2e  %4.2e \n'%(self.climate[i][0], self.climate[i][1], self.climate[i][2], self.climate[i][3], self.climate[i][4] , self.climate[i][5] ))
            fclimate.close()

        # Fallow Flag -   IFALLOW - Item 17
        ffallow = open(datapath + 'IFALLOW.in', 'w')
        for i in self.ifallow.keys():
            ffallow.write('%4d  %4d \n'%(i, self.ifallow[i] ))
        ffallow.close()

        # #
        # Crop Cost Parameters  -   [WPF-Slope, WPF-Int , Crop-Price] - Item 18
        fcost = open(datapath + 'CROPCOST.in', 'w')
        for i in self.cropcost.keys():
            fcost.write('%4d  %4f  %4f  %4f  \n'%(i, self.cropcost[i][0], self.cropcost[i][1], self.cropcost[i][2] ))
        fcost.close()

        ##
        # ALLOTGW    - for each stress period, and for each farm ALLOTGW  - Item 25
        if self.iallotgw == 2 and self.allotgw is not None:
            # for each SP
            for sp in self.allotgw.keys():
                fallotgw = open(datapath + 'ALLOTGW_%s.in'%str(sp), 'w')
                for fd in range(self.nfarms):
                    # [FarmID, GW_Allotment]
                    fallotgw.write('%4d  %4e \n'%(fd + 1, self.allotgw[sp][fd] ))
                fallotgw.close()

        # Non-routed deliveries file - for each stress period Item 36 (for each farm)
        if self.inrdfl == 1 and self.nonrouted is not None:
            # for each SP
            for sp in self.nonrouted.keys():
                fnrouted = open(datapath + 'NONROUTED_%s.in'%str(sp), 'w')
                for fd in range(self.nfarms):
                    # [FarmID, NRDV , NRDR, NRDU]
                    listnr = self.nonrouted[sp][fd]
                    fnrouted.write('%4d  %4e %4d %4d \n'%(fd + 1, listnr[0], listnr[1], listnr[2]  ))
                fnrouted.close()

        # ALLOTSW    - for each stress period, and for each farm ALLOTSW  - Item 39
        if self.iallotsw == 2 and self.allotsw is not None:
            # for each SP
            for sp in self.allotsw.keys():
                fallotsw = open(datapath + 'ALLOTSW_%s.in'%str(sp), 'w')
                for fd in range(self.nfarms):
                    # [FarmID, GW_Allotment]
                    fallotsw.write('%4d  %4e \n'%(fd + 1, self.allotsw[sp][fd] ))
                fallotsw.close()

    def write_file(self):
        """
        Write the package file.

        Returns
        -------
        None

        """

        # Run the method for writing external file *.in
        self.write_in_files()

        # -open file for writing
        f = open(self.fn_path, 'w')

        # -write header
        f.write('{}\n'.format(self.heading))

        # Item 1 - [PARAMETER NPFWL MXL {MXLP}]
        f.write('PARAMETER'+ '{:2d}{:2d}\n'.format(self.npfwl, self.mxl))
        # -Section 1 - FLAG BLOCKS  - flags are to be specified by blocks.
        f.write('FLAG_BLOCKS\n')

        # Item 2c -  MXACTW NFARMS NCROPS NSOILS
        f.write('{:4d}{:4d}{:4d}{:4d}\n'.format(self.mxactw, self.nfarms, self.ncrops, self.nsoils))

        # Item 2c - IFRMFL IRTFL ICUFL IPFL IFTEFL IIESWFL IEFFL
        f.write('{:4d}{:4d}{:4d}{:4d}{:4d}{:4d}{:4d}\n'.format(self.ifrmfl, self.irtfl, self.icufl, self.ipfl, self.iftefl, self.iieswfl, self.ieffl))

        # Item 2c  - IEBFL IROTFL IDEFFL IALLOTGW (WATER POLICY FLAGS)
        f.write('{:4d}{:4d}{:4d}{:4d}\n'.format(self.iebfl, self.irotfl, self.ideffl, self.iallotgw))

        # Item 2c - ICCFL
        f.write('{:4d}\n'.format(self.iccfl))

        # Item 2c - INRDFL {MXNRDT} ISRDFL IRDFL ISRRFL IRRFL IALLOTSW {PCLOSE}
        # if self.inrdfl == 0 and self.iallotsw <= 1:
        #     f.write('{:4d}{:4d}{:4d}{:4d}{:4d}{:4d}\n'.format(self.inrdfl, self.isrdfl, self.irdfl, self.isrrfl, self.irrfl, self.iallotsw))
        # if self.inrdfl == 1 and self.iallotsw <= 1:
        #     f.write('{:4d}{:4d}{:4d}{:4d}{:4d}{:4d}{:4d}\n'.format(self.inrdfl, self.mxnrdt,self.isrdfl, self.irdfl, self.isrrfl, self.irrfl, self.iallotsw))
        # if self.inrdfl == 0 and self.iallotsw > 1:
        #     f.write('{:4d}{:4d}{:4d}{:4d}{:4d}{:4d}{:4f}\n'.format(self.inrdfl, self.isrdfl, self.irdfl, self.isrrfl, self.irrfl, self.iallotsw,self.pclose))
        # if self.inrdfl == 1 and self.iallotsw > 1:
        #     f.write('{:4d}{:4d}{:4d}{:4d}{:4d}{:4d}{:4d}{:4e}\n'.format(self.inrdfl, self.mxnrdt,self.isrdfl, self.irdfl, self.isrrfl, self.irrfl, self.iallotsw,self.pclose))

        if self.inrdfl == 0 :
            f.write('{:4d}{:4d}{:4d}{:4d}{:4d}{:4d}'.format(self.inrdfl, self.isrdfl, self.irdfl, self.isrrfl, self.irrfl, self.iallotsw))
        if self.inrdfl == 1 :
            f.write('{:4d}{:4d}{:4d}{:4d}{:4d}{:4d}{:4d}'.format(self.inrdfl, self.mxnrdt,self.isrdfl, self.irdfl, self.isrrfl, self.irrfl, self.iallotsw))

        # IF iallotsw > 1 , PCLOSE must be specified:
        if self.iallotsw <= 1:
            f.write('\n')
        else:
            f.write('  %4f \n'%self.pclose)

        # Item 2c - IFWLCB IFNRCB ISDPFL IFBPFL IETPFL IRTPFL
        f.write('{:4d}{:4d}{:4d}{:4d}{:4d}{:4d}'.format(self.ifwlcb, self.ifnrcb, self.isdpfl, self.ifbpfl, self.ietpfl, self.irtpfl ))

        # Item2c (optional flags): {IRTPFL} {IOPFL} {IPAPFL}
        if self.isrdfl > 0:
            f.write('  %4d'%self.irtpfl)
        if self.ideffl > 0:
            f.write('  %4d'%self.iopfl)
        if self.iallotsw > 0:
            f.write('  %4d'%self.ipapfl)
        #-
        f.write('\n')


        # flag for auxiliary variables
        f.write('NOAUX\n')

        # flag for options
        f.write('NOOPT\n')

        # Items 3 to 4 if npfwl > 1
        # if npfwl > 1:
        #
        # Item 5  - GSURF
        f.write('OPEN/CLOSE DATA_FMP/GSE.in   1.0 (FREE) 1 \n')
        # Item 6 - FID
        if self.ifrmfl == 1:
            f.write('OPEN/CLOSE DATA_FMP/FID.in   1   (FREE) 1 \n')
        # Item 7 - OFE
        if self.ieffl == 1:
            f.write('OPEN/CLOSE DATA_FMP/OFE.in  \n')
        # Item 8 - SID
        f.write('OPEN/CLOSE DATA_FMP/SID.in   1 (FREE) -1 \n')
        # Item 9 - SOIL list
        f.write('OPEN/CLOSE DATA_FMP/SOILLIST.in  \n')
        # Item 10 - CID
        if self.irotfl >= 0:
            f.write('OPEN/CLOSE DATA_FMP/CID.IN   1 (FREE) -1 \n')
        # Item 11 - ROOT
        if self.irtfl == 1:
            f.write('OPEN/CLOSE DATA_FMP/ROOT.in  \n')
        # Item 12  - FTR FEP FEI
        if self.iftefl == 1:
            f.write('OPEN/CLOSE DATA_FMP/CROPEVAP.in \n')
        # Item 13 - FIESWP FIESWI
        if self.iieswfl == 1:
            f.write('OPEN/CLOSE DATA_FMP/CROPFIES.in \n')
        # Item 14 - PSI
        if self.iccfl == 1 or self.iccfl == 3 :
            f.write('OPEN/CLOSE DATA_FMP/PSI.in \n')
        # Item 15 - GROWTH
        if self.irtfl == 3 or self.icufl == 3 or self.ipfl == 3:
            f.write('OPEN/CLOSE DATA_FMP/GROWTH.in \n')
        # Item 16 - Climate Data
        if self.irtfl == 3 or self.icufl == 3 or self.ipfl == 3:
            f.write('OPEN/CLOSE DATA_FMP/CLIMATE.in \n')
        # Item 17 - IFALLOW
        if self.ideffl == -2 :
            f.write('OPEN/CLOSE DATA_FMP/IFALLOW.in \n')
        # Item 18 - CROP COST
        if self.ideffl > 0 and self.iben == 1 :
            f.write('OPEN/CLOSE DATA_FMP/CROPCOST.in \n')
        # Item 19 - FARM COST
        if self.ideffl > 0 and self.icost == 1 :
            f.write('OPEN/CLOSE DATA_FMP/FARMCOST.in \n')
        # Item 20 - ALLOTGW
        if self.iallotgw == 1 :
            f.write('OPEN/CLOSE DATA_FMP/ALLOTGW.in \n')
        # Item 21a - Semi-Routed SW
        if self.isrdfl == 1 :
            f.write('OPEN/CLOSE DATA_FMP/SRSWLOC.in \n')

        # Data for each stress period
        for p in range(self.nper):
            # To be correct : NP  = 0 set by default for each stress period
            f.write('{:10d}{:10d}\n'.format(self.itmp[p], 0))
            # Repeat this section ITMP times -  Dataset  23
            # Layer Row Column Farm-Well-ID Farm-ID QMAX
            if self.itmp[p] > 0:
                for well in range(self.itmp[p]):
                    wlist = self.farmwells[p][well]
                    f.write('{:4d}{:4d}{:4d}{:4d}{:4d}  {:5f}\n'.format(wlist[0]+1, wlist[1] +1, wlist[2] + 1, wlist[3], wlist[4], wlist[5]) )

            # Dataset 25 - ALLOTGW GW Allotments specified for each stress period.
            if self.iallotgw == 2:
                f.write('OPEN/CLOSE DATA_FMP/ALLOTGW_%s.in \n'%str(p))

            # dataset 29 - ROOT for each crop for each sp
            if self.irtfl == 2:
                #f.write('OPEN/CLOSE DATA_FMP/ROOT_.in \n')            #GDF (see issue #194)
                f.write('OPEN/CLOSE DATA_FMP/ROOT_%s.in \n'%str(p))    #GDF (see issue #194)

            if self.icufl == -1:
                # dataset 30a - KC for each crop for each sp
                #f.write('OPEN/CLOSE DATA_FMP/KC.in \n')               #GDF (see issue #194)
                f.write('OPEN/CLOSE DATA_FMP/KC_%s.in \n'%str(p))      #GDF (see issue #194)
                # dataset 30b - ET ref array for each sp
                f.write('OPEN/CLOSE DATA_FMP/ET_%s.in   1.0 (FREE) -1\n'%str(p))
                #f.write('OPEN/CLOSE DATA_FMP/ET_%s.in  \n'%str(p))

            if self.ipfl == 2:
                # dataset 33 - Precipitation PFLX, for each stress period, as CONSTANT (lumped parameter)
                f.write('CONSTANT %4.2e\n'%self.precip[p])
            # datset 36 , IF non-routed deliveries exist
            if self.inrdfl == 1:
                # Farm-ID (NRDV NRDR NRDU)
                f.write('OPEN/CLOSE DATA_FMP/NONROUTED_%s.in \n'%str(p))

            # Data Set 39, farm-specific water rights calls list required for each SP
            if self.iallotsw == 2:
                # read file where farmId and Water Calls are written
                f.write('OPEN/CLOSE DATA_FMP/ALLOTSW_%s.in \n'%str(p))


        f.close()
