# ******************************************************
## Revision "$LastChangedDate: 2015-05-21 09:31:08 +0200 (Thu, 21 May 2015) $"
## Date "$LastChangedRevision: 507 $"
## Author "$LastChangedBy: arthurbeusen $"
## URL "$HeadURL: https://pbl.sliksvn.com/globalnutrients/N_model/trunk/tools/n_groundwater.py $"
# ******************************************************
# Python modules
import os
import math
# Own modules
import ascraster
from error import *
import my_sys
import make_fraction
import aggregate
import litho_class
import sgd
import simpson


# Parameter to minimize the amount of output to the screen
max_lines = 3000
ldebug = False
def print_debug(val):
    if (ldebug):
        print val

def interpolate(x,xlist,ylist):
    '''
    Interpolate for given x and returns the interpolated value of y.

    If x is out of the boundary of xlist, then first or last value is returned.
    '''
    if (x < xlist[0]):
        return ylist[0] 
        #raise MyError ("Interpolation not possible for x: "+str(x)+" because list begins at: "+str(xlist[0]))
    elif (xlist[-1] < x):
        return ylist[-1]
        #raise MyError ("Interpolation not possible for x: "+str(x)+" because list ends at: "+str(xlist[-1]))
    else:
        # Interpolation
        xbegin = xlist[0]
        xend   = xlist[-1]
        ibegin = 0
        iend = len(xlist)-1
        for item in range(1,len(xlist)):
            if (x >= xlist[item]):
                xbegin = xlist[item]
                ibegin = item
            else:
                break
        for item in range(len(xlist)-2,0,-1):            
            if (x < xlist[item]):
                xend = xlist[item]
                iend = item
            else:
                break
        # Now interpolation
        if (ibegin == iend):
            return ylist[ibegin]
        else:
            return ylist[ibegin] + (x - xlist[ibegin])*(ylist[iend] - ylist[ibegin])/(xlist[iend] - xlist[ibegin]) 

class Cell():
    '''Derived parameter values for a cell'''
    def __init__(self,current_year,Tr,k,hist_years,hist_vals,D):
        self.current_year = current_year 
        self.Tr = Tr
        self.k = k
        self.hist_years = hist_years
        self.hist_vals = hist_vals
        self.D = D
        
    def ngrw(self,x):
        '''
        Calculate the historical groundwater load as an approximation of an integral

        Transport, delay and decay of a solute in groundwater from top aquifer to
        draining water courses or river.

        Decay with time during transport along a stream line in a homogeneous and

        isotropic aquifer, drained by parallel rivers or streams, is described 

        by a first order degradation proces:
        y(t)=y0(t).exp(-kt)                    (1)

        Assume inflow-concentration with time at the top of the aquifer, y0(t),

        can be described as a lineair function between known historical points. (2)
        Age of the groundwater (t) is described as

        t = -Tr.ln(x)                          (4)

        where  x=1-(z/D)                         (5)

        D  = maximum depth from top aquifer (m)

        Tr = mean time of residence in aquifer (year)
        t     z   x
        0     0   1
        inf   D   0
        Outflow concentration at depth x:

        y(x) = y0(-Tr.ln(x)).exp(k.Tr.ln(x))   (6)

        Integration is done by approximation by spliting up the aquifer in small pieces and calculating the integral.
        Integration of y(x) over 0 < x < 1 yields the average y
        Sum of Y is returned.
        Simpson's rule is applied here:
        S = (h/3) *(y0 + 4y1 + 2y2 + 4y3 ... + 2y(n-2) + 4y(n-1) +yn
        h = width of subinterval and n must be even.
        Now addaptive Simpson is used. 
        '''
        try:
            age = -self.Tr * math.log(1.0 - x/self.D)
            time = self.current_year - age
            y0t = interpolate(time,self.hist_years,self.hist_vals)
            yt = y0t * math.exp(-self.k * age)
        except (ValueError,OverflowError):
            y0t =0.0
            yt = 0.0
        return yt

    def denitr(self,x):
        '''

        Calculate the historical denitrification.
    
        Decay with time during transport along a stream line in a homogeneous and
        isotropic aquifer, drained by parallel rivers or streams, is described 
        by a first order degradation proces:

        y(t)=y0(t).exp(-kt)                    (1)
        Assume inflow-concentration with time at the top of the aquifer, y0(t),
        can be described as a lineair function between known historical points. (2)
        Age of the groundwater (t) is described as
        t = -Tr.ln(x)                          (4)

        where  x=1-(z/D)                         (5)
        D  = maximum depth from top aquifer (m)
        Tr = mean time of residence in aquifer (year)
        t     z   x
        0     0   1

        inf   D   0
        Outflow concentration at depth x:
        y(x) = y0(-Tr.ln(x)).exp(k.Tr.ln(x))   (6)

        Total denitrification at depth x:

        DenTot(x) = y0(-Tr.ln(x)).(1 - exp(k.Tr.ln(x)))   (7)
    
        Denitrification in the current year:
        Den(x) = y0(-Tr.ln(x)).[exp(k.{Tr.ln(x)+1}) - exp(k.Tr.ln(x)))]  
               = y0(-Tr.ln(x)).[exp(k) - 1].[exp(k.Tr.ln(x))]  (8)

        The assumption of equation 8 is that age of the groundwater is greater than 1 year.
    
        In the case the age of the groundwater is younger than 1 year, then eq (8) must be:
        Den(x) = y0(-Tr.ln(x)).[(1 - exp(k.Tr.ln(x)))]  (9)
    

        Groundwater of 1 year old is found on x1 = exp(-1/Tr)
    
        Eq (8) can be rewritten into:
        Den(x) = y0(-Tr.ln(x)).exp(k.Tr.ln(x))[exp(k) - 1]  x < x1 (10)
        Eq(10) says that denitritification in the current year is a constant factor (depending on first order decay rate) of the outflow 
        For the steady state situation, we can calculate the denitrification of the current year. Then y0 is constant in time.

    
        The integral for 0 <= x <= 1 is divided into 2 parts.
        For the part where the groundwater age is younger than 1 year:
    
        integral from x1 to 1 of eq(9) dx .

        Because the exp(k.Tr.ln(x)) = exp(ln(x))^(k.Tr) = x ^ (k.Tr)
        The primitive of function in eq (9) is y0.[x - 1/(k.Tr+1) * x^(k.Tr+1)]
        Fill in the boundary x=1 results in y0.[1-1/(kTr+1)] and x=x1 in y0.[x1 -  1/(k.Tr+1) * (x1)^(k.Tr+1)]
        So the result of part 1 of the integral is y0.[k.Tr/(kTr+1) - exp(-1/Tr) + 1/(k.Tr+1) * (exp(-1/Tr))^(k.Tr+1)]
    

        The second part of the integral is 
        integral from 0 to x1 of eq(10) dx = integral of (y0.(exp(k)-1) . x^(tr.k)
        The primitive of function in eq (10) is y0.(exp(k)-1).[ 1/(k.Tr+1) * x^(k.Tr+1)]
        Fill in the boundary x=0  results in 0 and for x=x1:  y0.(exp(k)-1).[1/(k.Tr+1) * (x1)^(k.Tr+1)]
        So the result of part 2 of the integral is y0.(exp(k)-1).[1/(k.Tr+1) * (x1)^(k.Tr+1)]

    
        Sum of the two parts gives the result of the whole integral (denitrification in the current year):
        Den = y0.k.Tr/(k.Tr+1) * [1 - exp(-1/Tr)] (11)
  
        Last simplification is because exp(k).x1^(Tr.k+1) = exp(k).exp(-1/Tr)^(Tr.k+1) = exp(k).exp((-Trk-1)/Tr) = exp(k).exp(-k-1/Tr) =

        exp(k-k-1/Tr) = exp(-1/Tr) 
     
        For the transient flow integration is done by approximation by spliting up the aquifer in small pieces and calculating the integral.  
        The denitrification function for this approximation uses a variation on equation (8):
        Den(x) = y0(-Tr.ln(x)).[exp(k.Tr.ln(x))].[exp(k.max{-Tr.ln(x),1}) - 1]  (12)

        Integration is done by approximation by spliting up the aquifer in small pieces and calculating the integral.
        Integration of y(x) over 0 < x < 1 yields the average y

        Sum of Y is returned.
        Simpson's rule is applied here:
        S = (h/3) *(y0 + 4y1 + 2y2 + 4y3 ... + 2y(n-2) + 4y(n-1) +yn
        h = width of subinterval and n must be even.
        Now addaptive Simpson is used. 
        '''
        try:
            age = -self.Tr * math.log(1.0 - x/self.D)
            time = self.current_year - age
            y0t = interpolate(time,self.hist_years,self.hist_vals)
            yt = y0t * (1.0 - math.exp(-self.k * age))
        except (ValueError,OverflowError):
            y0t =0.0
            yt = 0.0
        return yt

    def write(self):
        print self.current_year, self.Tr, self.k, self.hist_years, self.hist_vals, self.D

def calculate_steady_state(params,mask,\
              N_bal_grass,N_bal_arable,N_bal_natural,Qeff,fQgw,\
              total_balance,mouth_dict,basin):
    '''
    Groundwater transfer or delivery coefficient for Nitrate that is washed out
    from the soil and then is transported in a two-layered groundwater system to
    local water courses and distant rivers, depending on the depth and regional
    extent of the groundwater system. This program calculates mean residence or
    travel time, degradation rate coefficient and concentration in draining
    shallow and deep groundwater for the special case of steady state.
    Note: seepage flow is not taken into account
    Assume slow groundwaterflow in a regional deep aquifer of 50 m thickness and
    rapid subsurface interflow in shallow groundwater top-layer (thickness 5 m)
    Small water courses drain the shallow aquifer and large streams and rivers
    drain the deep aquifer. Both aquifers can be treated separately (parallel)
    or acting in a cascade of leaking (ground-)water reservoirs).
    '''
       
    if params.ldebug: print "Start with calculation of steady state shallow and deep groundwater concentrations."

    D_sgrw = params.D_sgrw  # [meter] Thickness of shallow groundwater layer
    D_dgrw = params.D_dgrw  # [meter] Thickness of deep groundwater layer

    # Read litho map and determine the effective porosity and dt50 for nitrogen
    litho_values = litho_class.Litho().make_value_list(params,mask,["porosity","dt50_shallow","dt50_deep"])

    # Read land area
    landarea = ascraster.Asciigrid(ascii_file=params.landarea,mask=mask,numtype=float)
    if params.ldebug: print params.landarea  + " has been read in calculation of steady state N groudwater." 

    # Make a copy (a deepcopy) of the return values
    N_deep_grw = ascraster.duplicategrid(Qeff)
    N_shallow_grw = ascraster.duplicategrid(Qeff)
    N_retention_grw = ascraster.duplicategrid(Qeff)
    N_top_dgrw = ascraster.duplicategrid(Qeff)
    N_denitr_grw = ascraster.duplicategrid(Qeff)
    N_denitr_grw_kgha = ascraster.duplicategrid(Qeff)
    
    for icell in range(landarea.length):
        area = landarea.get_data(icell,0.0)
        if (area > 0.0):
            N_bal = 0.0
            Qeff_cell = Qeff.get_data(icell,0.0)
            if (Qeff_cell > 0.0):
                # Get all the Nitrogen load
                N_bal += N_bal_grass.get_data(icell,0.0)
                N_bal += N_bal_arable.get_data(icell,0.0)
                N_bal += N_bal_natural.get_data(icell,0.0)
                
                if icell < max_lines: print_debug("N_bal: "+str(N_bal))
                # Convert load into application level (kg/km2).
                # N_bal = N_bal/area
                # Concentration of the inflow of the shallow groundwater (kg/1000 m3  = mg/l) because Qeff_cell is given in mm.
                # Cin_cell = N_bal/Qeff_cell
                poros = litho_values[icell][0]
                if icell < max_lines: print_debug("Porosity: "+str(poros))
                
                # Get partitioning of groundwater into shallow and deep.
                fQgw_cell = fQgw.get_data(icell,0.0)
                if icell < max_lines: print_debug("fQgw_cell: "+str(fQgw_cell))
                fQint_cell = 1.0 - fQgw_cell

         	# Mean time of residence of each groundwaterbody is calculated as follows:
        	# Tr=D*P/Qinflow ; D=thickness(m); P=eff. porosity(-); Qinflow=recharge(m/yr)
        	# vanDrecht et al. (2003), assumed a coupled system of a shallow and a deep
        	# groundwater reservoir, called a cascade of leaky groundwater reservoirs.                
                # The shallow groundwater partly recharges deep groundwater and discharges to
		# local shallow water courses and streams.
		# mean time of residence in the toplayer (eq. A15)
		# the factor 1000 is for converson of mm to m
	        # the limit of 1000 yr is for stability reasons
                # Residence time in shallow groundwater layer
                TRsgw=min(1000 * D_sgrw * poros/Qeff_cell,1000.0)
                if (fQgw_cell > 0.0):
                    TRdgw=min(1000 * D_dgrw * poros/(fQgw_cell * Qeff_cell),1000.0)
                else:
                    # No deep groundwater.
                    TRdgw = 1000.

                if icell < max_lines: print_debug("TRsgw: "+str(TRsgw))
                if icell < max_lines: print_debug("TRdgw: "+str(TRdgw))
                
		# Decomposition of Nitrate is assumed a first order process and the degradation
		# rate coefficient is calculated as follows (shallow groudwater is the same as deep groundwater):
		# Gerard did incorporate the effect of temperature on denitrification rate here by using the normalized
		# n_soil_denitrification.fr_denitrification_temperature. We don't do that here.
		DT50_shallow = litho_values[icell][1]
                DT50_deep    = litho_values[icell][2]
		if icell < max_lines: 
                    print_debug("DT50_shallow: "+str(DT50_shallow))
                    print_debug("DT50_deep: "+str(DT50_deep))
		
		if (DT50_shallow > 0.0):
		    Ksgw = math.log(2)/DT50_shallow
		    if icell < max_lines: print_debug("Ksgw calculated: "+str(Ksgw))
		else:
		    Ksgw = 0.0

		if (DT50_deep > 0.0):
		    Kdgw = math.log(2)/DT50_deep
		    if icell < max_lines: print_debug("Kdgw calculated: "+str(Kdgw))
		else:
		    Kdgw = 0.0

                if icell < max_lines: print_debug("Ksgw: "+str(Ksgw))
                if icell < max_lines: print_debug("Kdgw: "+str(Kdgw))

		# Special case is a steady state in the aquifer --> no concentration-change in
		# the inflow at the top boundary means agriculture is constant --> no trend in
		# fertilizer use (that is very special!)
		# The ratio f = outflow/inflow concentration of each groundwater layer can be
		# derived analytically for a homogeneous and mixed system:
		# f=1/(1+k.Tr), where k=degradation rate coefficient and Tr=time of residence
		# Nitrate outflow concentration is f times nitrate inflow concentration
		# at the top, and the fraction (1-f) is denitrificated.
		# delivery coefficient (DC) for shallow and deep groundwater
		DCsgw = 1./(1. + Ksgw * TRsgw)
		DCdgw = 1./(1. + Kdgw * TRdgw)
		if icell < max_lines: print_debug("DCsgw: "+str(DCsgw))
                if icell < max_lines: print_debug("DCdgw: "+str(DCdgw))

 
                # For practical implementation look at it is as follows:
                # Each kilo of Nitrate-N recharged to shallow groundwater is partitioned in:
                # flow(1) outflow to shallow groundwater draining water courses:   fQint*DCsgw
                # fQint=1-fQgw
                # flow(2) outflow to shallow groundwater draining deep groundwater: fQgw*DCsgw
                # flow(3) denitrification in shallow groundwater body: 1-DCsgw*(fQint+fQgw)
                RCsgw = 1. - DCsgw
                RCdgw = 1. - DCdgw

                # dgw flow(2) can be partitioned in:
                # 2a. outflow to deep groundwater draining rivers: fQgw*DCsgw*DCdgw
                # 2b. denitrification in deep groundwater: fQgw*DCsgw*RCdgw
                # mixed groundwater delivery coefficient:
                # outflow to draining water courses and rivers: 1+2a
                DCgw    = fQint_cell * DCsgw + fQgw_cell * DCsgw * DCdgw
                if icell < max_lines: print_debug("DCgw: "+str(DCgw))

                # Retention in shallow and deep groundwater together: 3+2b
                RCgw= RCsgw + fQgw_cell * DCsgw * RCdgw
                if icell < max_lines: print_debug("RCgw: "+str(RCgw))
                # where RC stands for Retention coefficient
                # Because the groundwater bodies are in steady state the kilos DCgw and
                # RCgw add up to 1

                # Calculate the N load of shallow groundwater and deep groundwater
                Nload_sgw      = N_bal * fQint_cell * DCsgw
                Nload_dgw      = N_bal * fQgw_cell * DCsgw * DCdgw
                Nload_top_dgrw = N_bal * fQgw_cell * DCsgw

                if icell < max_lines: print_debug("Nload_sgw: "+str(Nload_sgw))
                if icell < max_lines: print_debug("Nload_dgw: "+str(Nload_dgw))
                
                # Calculate the N denitrification of shallow groundwater and deep groundwater (total over whole period, steady state)
                N_retent_gw = N_bal * RCgw
                if icell < max_lines: print_debug("N_retent_gw: "+str(N_retent_gw))

                # Calculate the N denitrification of shallow groundwater emission of one year.
                N_denitr = N_bal * RCsgw
                N_denitr_kgha = N_denitr/(100.0 * area) # Conversion from km2 to ha (* 100.0) 

                # Store in the maps
                N_deep_grw.set_data(icell,Nload_dgw)
                N_shallow_grw.set_data(icell,Nload_sgw)
                N_retention_grw.set_data(icell,N_retent_gw)
                N_top_dgrw.set_data(icell,Nload_top_dgrw)
                N_denitr_grw.set_data(icell,N_denitr)
                N_denitr_grw_kgha.set_data(icell,N_denitr_kgha)

                if (math.fabs(N_bal - (Nload_dgw+Nload_sgw+N_retent_gw)) > 0.001):
                    print "ERROR: ", N_bal-(Nload_dgw+Nload_sgw+N_retent_gw),N_bal,Nload_dgw,Nload_sgw,N_retent_gw
            else:
                N_deep_grw.set_data(icell,0.0)
                N_shallow_grw.set_data(icell,0.0)
                N_retention_grw.set_data(icell,0.0)
                N_top_dgrw.set_data(icell,0.0)
                N_denitr_grw.set_data(icell,0.0)
                N_denitr_grw_kgha.set_data(icell,0.0)
        else:
            N_deep_grw.set_data(icell,0.0)
            N_shallow_grw.set_data(icell,0.0)
            N_retention_grw.set_data(icell,0.0)
            N_top_dgrw.set_data(icell,0.0)
            N_denitr_grw.set_data(icell,0.0)
            N_denitr_grw_kgha.set_data(icell,0.0)

    # Calculate the groundwater which flows directly to the coast (SGD)
    # Shallow and deep groundwater is reduced at the coast by SGD
    sgd.calculate(params,mask,mouth_dict,basin,N_shallow_grw,N_deep_grw)
        
    # Write rasters to output file
    N_deep_grw.write_ascii_file(os.path.join(params.outputdir,"N_deep_grw.asc"))
    N_shallow_grw.write_ascii_file(os.path.join(params.outputdir,"N_shallow_grw.asc"))
    N_retention_grw.write_ascii_file(os.path.join(params.outputdir,"N_retention_grw.asc"))
    # Write top of deep aquifer rasters to output file (can be used for the historical N directory, input of other years)
    N_top_dgrw.write_ascii_file(os.path.join(params.outputdir,"N_top_dgrw.asc"))
    # Write N denitrification of groundwater to output file
    N_denitr_grw.write_ascii_file(os.path.join(params.outputdir,"N_denitr_grw.asc"))
    N_denitr_grw_kgha.write_ascii_file(os.path.join(params.outputdir,"N_denitr_grw_kgha.asc"))

    # Aggregate rasters
    aggregate.aggregate_grid(basin,N_deep_grw,mouth_dict,item="ndgrw")
    aggregate.aggregate_grid(basin,N_shallow_grw,mouth_dict,item="nsgrw")
    aggregate.aggregate_grid(basin,N_retention_grw,mouth_dict,item="ngrw_retention")
    aggregate.aggregate_grid(basin,N_denitr_grw,mouth_dict,item="n_denitr_grw") 

    # Make groundwater loads as fraction of the total load
    temp = make_fraction.make_fraction(params,N_deep_grw,total_balance)
    # Write to output file
    temp.write_ascii_file(os.path.join(params.outputdir,"fr_N_deep_grw.asc"))
    aggregate.aggregate_grid(basin,temp,mouth_dict,item="fr_ndgrw",method=4,weighing=landarea)

    temp = make_fraction.make_fraction(params,N_shallow_grw,total_balance)
    # Write to output file
    temp.write_ascii_file(os.path.join(params.outputdir,"fr_N_shallow_grw.asc"))
    aggregate.aggregate_grid(basin,temp,mouth_dict,item="fr_nsgrw",method=4,weighing=landarea)    
    
    temp = make_fraction.make_fraction(params,N_retention_grw,total_balance)
    # Write to output file
    temp.write_ascii_file(os.path.join(params.outputdir,"fr_N_retention_grw.asc"))
    aggregate.aggregate_grid(basin,temp,mouth_dict,item="fr_ngrw_retention",method=4,weighing=landarea)
    
    # Delete the data which is not used anymore
    del landarea
  
    return N_deep_grw,N_shallow_grw,N_retention_grw                                       
    
def calculate_transient_state(params,mask,\
                              Qeff, fQgw,\
                              total_balance,mouth_dict,basin):
    '''
    Calculate the (historical) N load of groundwater
    The N_load of groundwater is returned as raster.
    This routine uses historical N leaching raster files,which are available on the directory params.N_histdir.

    Groundwater transfer or delivery coefficient for Nitrate that is washed out
    from the soil and then is transported in a two-layered groundwater system to
    local water courses and distant rivers, depending on the depth and regional
    extent of the groundwater system. This program calculates mean residence or
    travel time, degradation rate coefficient and concentration in draining
    shallow and deep groundwater for the special case of steady state.
    Note: seepage flow is not taken into account
    Assume slow groundwaterflow in a regional deep aquifer of 50 m thickness and
    rapid subsurface interflow in shallow groundwater top-layer (thickness 5 m)
    Small water courses drain the shallow aquifer and large streams and rivers
    drain the deep aquifer. Both aquifers can be treated separately (parallel)
    or acting in a cascade of leaking (ground-)water reservoirs).
    '''
    
    if params.ldebug: print "Start with calculation of transient shallow and deep groundwater concentrations."

    D_sgrw = params.D_sgrw  # [meter] Thickness of shallow groundwater layer
    D_dgrw = params.D_dgrw  # [meter] Thickness of deep groundwater layer

    # Read land area
    landarea = ascraster.Asciigrid(ascii_file=params.landarea,mask=mask,numtype=float)
    if params.ldebug: print params.landarea  + " has been read in calculation of transient state N groudwater." 

    # Read litho map and determine the effective porosity and dt50 for nitrogen
    litho_values = litho_class.Litho().make_value_list(params,mask,["porosity","dt50_shallow","dt50_deep"])
    
    # Fill in values for shallow groundwater and deep groundwater
    litho_values_shallow = []
    litho_values_deep = []
    for item in range(len(litho_values)):
        # Experiment with DT50 one year longer.
        #if (litho_values[item][1] > 0.0):
        #    litho_values_shallow.append([litho_values[item][0],litho_values[item][1]+1.0])
        #else:
        #    litho_values_shallow.append([litho_values[item][0],litho_values[item][1]])
        litho_values_shallow.append([litho_values[item][0],litho_values[item][1]])
        litho_values_deep.append([litho_values[item][0],litho_values[item][2]])

    # List with the historical years
    hist_years = []
    
    # Search for the historical directory.
    files = os.listdir(params.N_histdir)
    # Look for the year in the filename
    for file in files:
        for i in range(1500,2500):
            if (file.find(str(i)) > -1):
                # Found year
                lfound = False
                for y in range(len(hist_years)):
                    if (hist_years[y] == i):
                        lfound = True     
                        break
                if (not lfound):
                    hist_years.append(i)
                
                # Stop searching
                break
    # Take the current year as well
    lfound = False
    for y in range(len(hist_years)):
        if (hist_years[y] == params.year):
            lfound = True     
            break
    if (not lfound):
        hist_years.append(params.year)

    # Sort list of years:
    hist_years.sort()

    # Make a list of filenames for historical Nleaching files.
    nleaching_files = len(hist_years) * [""]
    top_dgrw_files = len(hist_years) * [""]

    # Try to put the filename into the right place in the array
    for file in files:
        iplace = -1
        for i in range(len(hist_years)):
            if (file.find(str(hist_years[i])) > -1):
                iplace = i
                print str(hist_years[i])," ",iplace," ",file
                break
        # Check whether year is found.
        if (iplace == -1):
            print "Historical file: " + os.path.join(params.N_histdir,file) + " is not used (no year info)!."
            # Go to next file
            continue
        # Check for leaching or shallow in filename
        if (file.upper().find("LEACH") > -1):
            # leaching is found
            nleaching_files[iplace] = os.path.join(params.N_histdir,file)
        elif (file.upper().find("TOP") > -1):
            # Top is found
            top_dgrw_files[iplace] = os.path.join(params.N_histdir,file)
        else:    
            print "Historical file: " + os.path.join(params.N_histdir,file) + " is not used (no category info)!."
            # Go to next file
            continue

    # Put leaching file of current year also in the array.
    iplace = -1
    for i in range(len(hist_years)):
        if (hist_years[i] == params.year):
            iplace = i
    nleaching_files[iplace] = os.path.join(params.outputdir,"Nleaching.asc")
    top_dgrw_files[iplace]  = os.path.join(params.outputdir,"N_top_dgrw.asc")
    
    # Do the historical calculation for shallow groundwater.
    N_shallow_grw,N_denitr_grw = do_N_history(params,mask,hist_years,nleaching_files,Qeff,D_sgrw,litho_values_shallow)
    
    # Aggregate raster
    aggregate.aggregate_grid(basin,N_denitr_grw,mouth_dict,item="n_denitr_grw")
    # Write raster to output file
    N_denitr_grw.write_ascii_file(os.path.join(params.outputdir,"N_denitr_grw.asc")) 

    # Calculate the N denitrification in kg/ha
    N_denitri_kgha = N_denitr_grw  # Make a grid for the output.
    for icell in range(N_denitr_grw.length):
        Nden_cell = N_denitr_grw.get_data(icell,0.0)
        area_cell = 100.0 * landarea.get_data(icell,0.0) # Conversion from km2 to ha
        try:
            Nden_cell_kgha = Nden_cell/area_cell
        except ZeroDivisionError:
            Nden_cell_kgha = 0.0
        N_denitri_kgha.set_data(icell,Nden_cell_kgha)

    N_denitri_kgha.write_ascii_file(os.path.join(params.outputdir,"N_denitr_grw_kgha.asc"))  

    # Calculate the N load at the top of the deep groundwater aquifer and outflow of the shallow groundwater.
    N_top_dgrw = ascraster.duplicategrid(N_shallow_grw)
    for icell in range(N_shallow_grw.length):
        fQgw_cell = fQgw.get_data(icell,0.0)
        load =  N_shallow_grw.get_data(icell,0.0)
        N_top_dgrw.set_data(icell,fQgw_cell * load)
        N_shallow_grw.set_data(icell,(1.0 - fQgw_cell) * load)
    
    # Write top of deep aquifer rasters to output file (can be used for the historical N directory, input of other years)
    N_top_dgrw.write_ascii_file(os.path.join(params.outputdir,"N_top_dgrw.asc"))
    
    # Remove this grid. Is not needed anymore.
    del N_top_dgrw
    
    # Calculate the deep groundwater flow
    # Loop over all the grid cells.
    water_top_dgrw = ascraster.duplicategrid(Qeff)
    for icell in range(Qeff.length):
        Qeff_cell = Qeff.get_data(icell,0.0)
        fQgw_cell = fQgw.get_data(icell,0.0)
        water_top_dgrw.set_data(icell,Qeff_cell*fQgw_cell)

    # Do the historical calculation for deep groundwater. Denitrification is calculated in do_N_history, but for deep groundwater, we do not use this.
    N_deep_grw,N_denitr_grw = do_N_history(params,mask,hist_years,top_dgrw_files,water_top_dgrw,D_dgrw,litho_values_deep)
    
    # Calculate the groundwater which flows directly to the coast (SGD)
    # Shallow and deep groundwater is reduced at the coast by SGD
    sgd.calculate(params,mask,mouth_dict,basin,N_shallow_grw,N_deep_grw)

    # Write rasters to output file
    N_shallow_grw.write_ascii_file(os.path.join(params.outputdir,"N_shallow_grw.asc"))
    N_deep_grw.write_ascii_file(os.path.join(params.outputdir,"N_deep_grw.asc"))

    # Aggregate rasters
    aggregate.aggregate_grid(basin,N_shallow_grw,mouth_dict,item="nsgrw")
    aggregate.aggregate_grid(basin,N_deep_grw,mouth_dict,item="ndgrw")

    # Make groundwater loads as fraction of the total load and write to output file
    temp = make_fraction.make_fraction(params,N_shallow_grw,total_balance)
    temp.write_ascii_file(os.path.join(params.outputdir,"fr_N_shallow_grw.asc"))
    aggregate.aggregate_grid(basin,temp,mouth_dict,item="fr_nsgrw",method=4,weighing=landarea)    

    temp = make_fraction.make_fraction(params,N_deep_grw,total_balance)
    temp.write_ascii_file(os.path.join(params.outputdir,"fr_N_deep_grw.asc"))
    aggregate.aggregate_grid(basin,temp,mouth_dict,item="fr_ndgrw",method=4,weighing=landarea)

    # To be consistent with the steady state approach, we "calculate"  the retention in groundwater.
    # Because of the transient historical flow, we can not compute this.
    # So we make everything zero.

    N_retention_grw = ascraster.duplicategrid(N_deep_grw)
    N_retention_grw.add_values(N_retention_grw.length * [0.0])

    # Aggregate rasters
    aggregate.aggregate_grid(basin,N_retention_grw,mouth_dict,item="ngrw_retention")

    # Make groundwater retention as fraction of the total load and write to output file   
    temp = make_fraction.make_fraction(params,N_retention_grw,total_balance)
    temp.write_ascii_file(os.path.join(params.outputdir,"fr_N_retention_grw.asc"))
    aggregate.aggregate_grid(basin,temp,mouth_dict,item="fr_ngrw_retention",method=4,weighing=landarea)
    
    return N_deep_grw,N_shallow_grw,N_retention_grw
    
            
def do_N_history(params,mask,hist_years,hist_files,Qeff,D,litho_values):
    '''
    Do the calculation of the historical N groundwater load for each grid cell.
    '''
    
    ldebug = params.ldebug
    
    if (ldebug == 1): print "Start of do_N_history."
    # The maximum of 50 years is taken for the history of N. This number is taken because with a DT50 for nitrogen of 2 years,
    # the decay rate is ln(2)/dt50 = 0.3466, the reduction of exp(-kt) for t-50  is almost 3e-8.
    # The nitrogen load of more than 50 years old is assumed to be zero. Beautifull story, but we use 150 years!
    
    #if (ldebug == 1): print "A maximum of 50 years is taken for the history."
    if (ldebug == 1): print "Current year is " + str(params.year)
    
    if (ldebug == 1): 
        print "Years and file names of historical N leaching files."
        for i in range(len(hist_years)):
            print hist_years[i], hist_files[i]
    
    current_year = params.year        
    
    # Find the current year
    iyear_end = len(hist_years)-1
    iyear_begin  = 0
    for item in range(1,len(hist_years)):
        # Take all the historical files which are at most 150 years old.
        if (current_year-150 >= hist_years[item]):
            iyear_begin = item
        else:
            break
    for item in range(1,len(hist_years)):        
        for item in range(len(hist_years)-2,0,-1): 
            if (current_year <= hist_years[item]):
                iyear_end = item
            else:
                break
    if (ldebug == 1):  print "Start periode for do_N_history is: ",hist_years[iyear_begin]," and ",hist_years[iyear_end]
    # Start with the initialisation
    year_end = hist_years[iyear_end]
    year_begin = hist_years[iyear_begin]

    # Read all the historical N load files for groundwater.
    hist_grid = []
    years = []
    counter = 0
    for iyear in range(iyear_begin,iyear_end+1):
        hist_grid.append(ascraster.Asciigrid(ascii_file=hist_files[iyear],mask=mask,numtype=float))
        if (ldebug == 1):
            print str(hist_files[iyear]) + " is read in calculation of do_N_history."
    
        # Compare whether the grid has the right dimension
        equal = hist_grid[0].compare(hist_files[iyear])
        if (not equal):
            raise MyError("Historical N leaching files don't have the same header. Files: " + \
                           hist_files[iyear] + " and " + hist_files[iyear_begin])
                                            
        years.append(hist_years[iyear])
        counter += 1
    
    # Make a copy (a deepcopy) of the return values
    N_grw = ascraster.duplicategrid(Qeff)
    N_denitr_grw = ascraster.duplicategrid(Qeff)

    import time
    start_time = time.time()
    try:
        import pp
        lpp  = True
        if (params.ncpu == 1):
            # No parallelization of the groundwater calculation
            lpp = False
            print "lpp is FALSE because ncpu == 1"
        else:
            print "lpp is TRUE"
    except ImportError:
        lpp = False
        print "lpp is FALSE"

    if (lpp):
        if (params.ncpu == 1):
            # No parallelization of the groundwater calculation
            lpp = False

    if (lpp):
        print "PP is available, so parallel computing is used for groundwater."
        # Make preparations for parallel computing.
        # Initialisation of local host.
        #job_server = pp.Server(ncpus=24,ppservers=("131.224.176.58","131.224.176.29"))
        #job_server = pp.Server(ppservers=("*:54321","localhost"),secret="2ICT4PBL")
        job_server = pp.Server()
        
        # Initialisation of the number of CPUs
        if (params.ncpu < 0):
            # Take all the cores available in the computer
            job_server.set_ncpus()
        else:
            if (job_server.get_ncpus() < params.ncpu):
                job_server.set_ncpus()
            else:    
                job_server.set_ncpus(params.ncpu)
        print "PARALLEL: There are ",job_server.get_ncpus()," cores available."

        # Initalisation of job array
        jobs = []
        max_parts = 63
        start = 0
        end = N_grw.length
        step = max(N_grw.length / max_parts,500)
        parts = 1 + int(end/step) 
        grw_all = []
        denit_all = []
        for index in range(parts):
            starti = start+index*step
            endi = min(start+(index+1)*step, end)
            print " Submit block ",index," ",starti,"  ",endi
            #        (math.sin, starti,endi),(),("math",)))
            jobs.append(job_server.submit(grw_cell,\
                     (starti,endi,Qeff,litho_values,\
                     D,counter,hist_grid,current_year,years,params.epsilon_grw),\
                     (Cell,simpson.adaptive_simpsons_rule,simpson.recursive_asr,\
                      simpson.simpsons_rule,interpolate),("simpson","math")))

        # Get all results from all jobs.    
        for job in jobs:
            grw,denit = job()
            grw_all.extend(grw)
            denit_all.extend(denit)
        # Print job statistics to screen.
        job_server.print_stats()
        #job_server.destroy()

    else:
        print "PP is NOT available, so serial computing (one core) is used for groundwater."
        # Calculate groundwater outflow and groundwater denitrification for each grid cell.
        grw_all = []
        denit_all = []
        icell_start = 0
        icell_end = Qeff.length
        grw_all,denit_all = grw_cell(icell_start,icell_end,Qeff,litho_values,\
                                     D,counter,hist_grid,current_year,years,params.epsilon_grw)

    end_time = time.time()
    print "Wall clock time groundwater: ",end_time - start_time," seconds."
    N_grw.add_values(grw_all) 
    N_denitr_grw.add_values(denit_all) 

    return N_grw, N_denitr_grw


def grw_cell(icell_start,icell_end,Qeff,litho_values,D,counter,hist_grid,current_year,years,eps):
    '''
    Calculation of transient groundwater for one cell.
    Result of this function is put in lists grw and denit
    '''
    grw = []
    denit = []
    for icell in xrange(icell_start,icell_end):
        Qeff_cell = Qeff.get_data(icell,0.0)
        if (Qeff_cell > 0.0):
            poros = litho_values[icell][0]                    
            # Mean time of residence of each groundwaterbody is calculated as follows:
            # Tr=D*P/Qinflow ; D=thickness(m); P=eff. porosity(-); Qinflow=recharge(m/yr)
            # vanDrecht et al. (2003), assumed a coupled system of a shallow and a deep
            # groundwater reservoir, called a cascade of leaky groundwater reservoirs.                
            # The shallow groundwater partly recharges deep groundwater and discharges to
            # local shallow water courses and streams.
            # mean time of residence in the toplayer (eq. A15)
            # the factor 1000 is for converson of mm to m
            # the limit of 1000 yr is for stability reasons
            # Residence time in shallow groundwater layer
            Tr =min(1000 * D * poros/Qeff_cell,1000.0)
                      
            # Decomposition of Nitrate is assumed a first order process and the degradation
            # rate coefficient is calculated as follows (shallow groudwater is the same as deep groundwater):
            DT50 = litho_values[icell][1]
            if (DT50 > 0.0):
                # Decomposition of Nitrate is active
                K = math.log(2)/DT50
            else:
                # No decomposition of Nitrate
                K = 0.0
            # Get all historical load for this cell
            Nload_top_layer = []
            for i in range(counter):
                val = hist_grid[i].get_data(icell,0.0)
                Nload_top_layer.append(val)
	            
            # Calculation of the historical groundwater load
            cell = Cell(current_year,Tr,K,years,Nload_top_layer,D)            
            if (Tr > 0.0):    
                N_grw_cell = simpson.adaptive_simpsons_rule(cell.ngrw,0.0,D,eps)
                # Calculate the average load of the aquifer.       
                N_grw_cell /= D
            else:
                # Tr = 0, so also no decay
                N_grw_cell = interpolate(current_year,years,Nload_top_layer)

            # Store in the map
            #N_grw.set_data(icell,N_grw_cell)               

            # Calculation of the historical denitrification of shallow groundwater
            if (Tr > 0.0):               
                N_denitr_grw_cell = simpson.adaptive_simpsons_rule(cell.denitr,0.0,D,eps)
                # Calculate the average load of the aquifer.       
                N_denitr_grw_cell /= D
            else:
                # Tr = 0, so also no decay
                N_denitr_grw_cell = 0.0

            # Store denitrifcation in the map.
            #N_denitr_grw.set_data(icell,N_denitr_grw_cell)  

            # Check steady state versus transient (only for checking purpose in the year 1900!)
            #fr_outflow_ss = 1.0/(1.+Tr*K)
            #N_grw_ss = interpolate(current_year,years,Nload_top_layer) * fr_outflow_ss
            #N_denitr_grw_ss = interpolate(current_year,years,Nload_top_layer) * (1.0 - fr_outflow_ss)
            #if (math.fabs(N_grw_ss - N_grw_cell) > 10.0):
            #    print "ERROR TRSGRW: ", icell,N_grw_ss, N_grw_cell,N_grw_ss - N_grw_cell, Tr,K, Nload_top_layer, N_denitr_grw_ss, N_denitr_grw_cell,cell.write()
            #if (math.fabs(N_denitr_grw_ss - N_denitr_grw_cell) > 10.0):
            #    print "ERROR DENGRW: ", icell,N_denitr_grw_ss, N_denitr_grw_cell,N_denitr_grw_ss - N_denitr_grw_cell, Tr,K, Nload_top_layer, N_grw_ss, N_grw_cell,cell.write()
	    	    
        else:
            N_grw_cell = 0.0
            N_denitr_grw_cell = 0.0
        
        # Put result in output arrays
        grw.append(N_grw_cell)
        denit.append(N_denitr_grw_cell)
    return grw, denit
