# ******************************************************
## 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_soil_denitrification.py $"
# ******************************************************

import os
import ascraster
import math
import make_fraction
import aggregate
from error import *

normalize_denitri_10_degree = 0.1511464474
ldebug = False

#oc_xlist = [0.0,1.0,3.0,6.0,100.0]
#oc_ylist = [0.0,0.1,0.2,0.3,  0.3]

oc_xlist = [0.0,6.0,100.0]
oc_ylist = [0.0,0.3,  0.3]

#ph_xlist=[0.0,3.0,4.0 ,4.5 ,5.0,5.5 ,6.0 ,6.5  ,7.0  ,7.5  ,8.0 ,8.5 ,9.0 ,9.5 ,10.0,11.0,100.0]
#ph_ylist=[1.0,1.0,0.83,0.77,0.5,0.35,0.20,0.087,0.045,0.087,0.20,0.35,0.50,0.77,0.83, 1.0,1.0  ]
ph_xlist=[0.0,3.0,4.0 ,4.5 ,5.0,5.5 ,6.0 ,6.5  ,7.0  ,7.5  ,100.0]
#ph_ylist=[1.0,1.0,0.83,0.77,0.5,0.35,0.20,0.087,0.045,0.0,0.0  ]
ph_ylist=[1.0,1.0,0.83,0.77,0.4,0.25,0.10,0.037,0.0,0.0,0.0  ]

ph1_xlist=[0.0,3.0,4.0 ,5.0 ,6.0 ,7.0  ,7.5  ,100.0]
ph1_ylist=[0.0,0.0,0.1 ,0.25,0.80,1.0,1.0,1.0  ]



def interpolate(x,xlist,ylist):
    '''
    Interpolate for given x and returns the interpolated value of y.
    '''
    if (x < xlist[0]):
        raise MyError ("Interpolation not possible for x: "+str(x)+" because list begins at: "+str(xlist[0]))
    elif (xlist[-1] < x):
        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]) 

def fr_Ndenitrification_ph(ph):
    '''
    Reduction of N denitrification based on pH value of the soil.
    The function is used in the riparian zone as a multiplier for the N2O emission.
    '''
    if (ph < 0.0): return 1.0
    return interpolate(ph,ph_xlist,ph_ylist)


def fr_Ndenitrification_temperature(temp):
    '''
    For the soil denitrification temperature stress factor we use Arrhenius
    see the NLEAP-model ref. Shaffer, M.J., Halvorson, A.D., Pierce, F.J. in:
    Managing Nitrogen for Groundwater Quality and Farm Profitability. Proc.
    symp. Crop Science Society of America, Anaheim, CA, 1988, eds. Follet,
    R.F., Keeney, D.R., Cruse, R.M.
    in SI: TFAC = A.EXP(-Ea/(R.K))
    A  = 1.68E9 constant
    Ea = activation energy = 13.0 Kcal/mol (=13.0*4.1868E3 J/mol = 54428.4)
    R  = universal gasconstant 8.3144 J/(K.MOL)
    K  = temperature in Kelvin (Celcius+273)
    calculate denitrification rate limited by temperature
    f_denitrification = 1.68e9*exp(-54428/(8.3144*(temperature+273)))

    Denitrification rate is normalized with a temperature of 10 degrees Celcius (factor 0.1511464474).
    Why this is normalized, i don't know. So i leave that out or not?!  
    '''
    #return (1.68e9*math.exp(-54428/(8.3144*(temp+273))))/normalize_denitri_10_degree
    # Temperature experiment with 1 degree higher.
    #return 7.94e12*math.exp(-74830/(8.3144*(temp+274)))
    return min(7.94e12*math.exp(-74830/(8.3144*(temp+273))),1.0)

def fr_Ndenitrification_oc(oc):
    '''
    fr_denitr_oc = exp(-0.094*1.72*oc)
    Note: 0.58 = OC/OM of C/H; in NLES wordt humus (organic matter) gebruikt.
    OM = 1.72 OC, maar dat geldt voor stabiele humus.
    in Van Drecht et al. werd ook nog het effect van drainage van het
    bodemprofiel op denitrificatie meegenomen, maar dat lijkt voldoende
    afgedekt te zijn met het effect van de organische stof.
    '''

    if (oc < 0.0): return 0.0
    return interpolate(oc,oc_xlist,oc_ylist)
    #if (oc < 1.0):
    #    return 0.0
    #elif (oc < 3.0):
    #    return 0.1
    #elif (oc < 6.0):
    #    return 0.2
    #else:
    #    return 0.3

def fr_Ndenitrification_text(text):
    '''
    fr_denitr_clay = exp(-0.03*clay)
    No explanation is given by van Drecht.
    '''
    if (text < 1):
        return 0.0
    elif (text < 1.5):
        # Coarse texture
        return 0.0
    elif (text < 2.5):
        # Medium texture
        return 0.1
    elif (text < 3.5):
        # Fine texture
        return 0.2
    elif (text < 4.5):
        # Very fine texture
        return 0.3
    elif (text < 5.5):
        # Organic texture (peat)
        return 0.0
    else:
        # Not known texture class
        return 0.0

def fr_Ndenitrification_drainage(drain):
    '''
    '''
    if (drain < 1):
        return 0.0
    elif (drain < 1.5):
        # Very poorly drained
        return 0.4
    elif (drain < 2.5):
        # poorly drained
        return 0.3
    elif (drain < 3.5):
        # Imperfectly drained
        return 0.2
    elif (drain < 4.5):
        # Moderately well drained
        return 0.1
    elif (drain < 7.5):
        # All other well drained
        return 0.0
    else:
        # Not known drainage class
        return 0.0


def trzd(Qeff,tawcd):
    '''
    Trz kan dus erg klein worden, maar dat geeft niet.
    In droge gebieden kan het neerslagoverschot zeer klein worden en daarmee
    de verblijftijd, Trz, erg groot. Gelukkig maakt het weinig uit of Trz
    groot of zeer groot is. In beide gevallen is uitspoeling verwaarloosbaar.
    Default 100 years when Qeff is very small.
    '''
    if (Qeff > 0.0):
        return tawcd/Qeff
    else:    
        return 100.0

    
def calculate_N_soil_denitrification(params,mask,grassarea,\
                       croparea,natarea,\
                       N_bal_grass,N_bal_arable,N_bal_natural,\
                       fQsro,pnet,total_balance,\
                       mouth_dict,basin):
    '''
    Calculates the N soil leaching and the N soil dentrification in the top soil (1 meter).
    It returns the ascraster object with the N soil dentrification.
    N_leaching is substracted from N_bal_grass, N_bal_arable and N_bal_natural.
    
    Soil N leaching model of 2004, based on Drecht 2002 and the Danish NLES-
    model (ref. see Euroharp website)
    Leaching of nitrate from the soil compartment to shallow groundwater is
    calculated from the extended surface N balance:
    surface balance N surplus, N surface runoff and soil N denitrification:
    Nsurp = Ninputs - Noutputs(including N surface runoff)
    Nlea = Nsurp - Nden = Nsurp - Fden*Nsurp [kg/ha/yr],
    Nlea/Nsurp = Fleach = 1 - Fden
    where Fden is the fraction of the Nsurp that is denitrificated.
    Denitrification of nitrate-N depends on the degree of aeration,
    that is related to the water content and the texture: fine texture means
    more water means less fresh air means that denitrification is favoured.
    High temperatures favors mineralization and denitrification while
    low temperatures stop all N transformation processes.
    In dry climates flow velocities in the soil compartment are small, so
    nitrate in the soil moisture will stay in the soil compartment for a
    long (residence)time. In humid and wet conditions the risk of leaching
    of the mobile nitrate is much higher.
    By definition the annual refresh rate of the mobile soil moisture in
    the soil compartment is the number of mobile moisture volumes (MV)
    that is replaced by fresh rainfall excess:
    MV = rainfall_excess/tawc [1/yr],
    which is the reciproce value of the mean residence time or travel time
    of water in the soil, Trz [yr].
    tawc = Total Available soilWater Capacity in the rooting zone(mm)
    in fact we have tawc-maps for dominant, suitable and natural landuse,
    in this case we use a generic tawc for suitable soils
    pnet is voor rainfed landbouw en natuur (notcultivated)
    maar voor cropland moeten we rekening houden met irrigatie (ariede gebieden)
    en kunnen we dus beter Qeff (zie WBAL).
    
    Nlea/Nsurp = Fuse*exp(-a*Ftrz*Ftmp-Foc-Fclay); a=calibration factor
    
    Foc and Fclay are additional factors for organic carbon and clay content.
    Zeer kwetsbaar voor nitraatuitspoeling zijn arable land of cropland,
    no permanent crops, no under- or after crop, die een kort groeiseizoen
    hebben en daardoor ook een beperkt wortelstelsel. Mais is een rijenteelt
    en aardappelen worden geteelt op ruggen. De afstanden tussen de rijen en 
    ruggen is wel 60 cm. De zaak is dus maar half bedekt gedurende een deel
    van het groeiseizoen, dat pas laat in het voorjaar begint.
    de gewascoefficienten zijn ontleend aan het NLES model.
    '''
       
    if params.ldebug: print "Start with calculation of fraction soil leaching."

    # Gerard assumes that there is a DT50 of 2 years, therefore the
    # calibration factor a is equal to log(2)/2 = 0.3466
    a = 0.5 * math.log(2)

    # Read air temperature map
    temp = ascraster.Asciigrid(ascii_file=params.temperature,mask=mask,numtype=float)
    if params.ldebug: print "Air temperature is read in calculation of fraction soil leaching."
    
    tawcd = ascraster.Asciigrid(ascii_file=params.tawcd,mask=mask,numtype=float)
    if params.ldebug: print "Tawcd is read in calculation of fraction soil leaching."
       
    tawcs = ascraster.Asciigrid(ascii_file=params.tawcs,mask=mask,numtype=float)
    if params.ldebug: print "Tawcs is read in calculation of fraction soil leaching."

    s_oc = ascraster.Asciigrid(ascii_file=params.s_oc,mask=mask,numtype=float)
    if params.ldebug: print "S_Oc is read in calculation of fraction soil leaching."
    
    d_oc = ascraster.Asciigrid(ascii_file=params.d_oc,mask=mask,numtype=float)
    if params.ldebug: print "D_Oc is read in calculation of fraction soil leaching."

    s_txt = ascraster.Asciigrid(ascii_file=params.s_txt,mask=mask,numtype=int)
    if params.ldebug: print "S_texture is read in calculation of fraction soil leaching."
    
    d_txt = ascraster.Asciigrid(ascii_file=params.d_txt,mask=mask,numtype=int)
    if params.ldebug: print "D_texture is read in calculation of fraction soil leaching."

    s_drain = ascraster.Asciigrid(ascii_file=params.s_drain,mask=mask,numtype=int)
    if params.ldebug: print "S_drainage is read in calculation of fraction soil leaching."
    
    d_drain = ascraster.Asciigrid(ascii_file=params.d_drain,mask=mask,numtype=int)
    if params.ldebug: print "D_drainage is read in calculation of fraction soil leaching."


    # Make a copy (alias not a deepcopy) of the return values
    N_denitrification = temp
    Qeff = ascraster.duplicategrid(temp)
    fr_denitri = s_oc
    Nleach = s_txt
    if (Nleach.nodata_value == None):
        Nleach.nodata_value = -1
    if (fr_denitri.nodata_value == None):
        fr_denitri.nodata_value = -1
    # Make copy for agricultural and natural soil denitrification
    N_denitrification_arable  = ascraster.duplicategrid(temp)
    N_denitrification_grass   = ascraster.duplicategrid(temp)
    N_denitrification_nat     = ascraster.duplicategrid(temp)

    if (ldebug):
        fr_climate = ascraster.duplicategrid(temp)
        fr_oc = ascraster.duplicategrid(temp)
        fr_txt = ascraster.duplicategrid(temp)
        fr_drain = ascraster.duplicategrid(temp)
        Ftrzd = ascraster.duplicategrid(temp)
        Ftrzs = ascraster.duplicategrid(temp)
        Ftemp = ascraster.duplicategrid(temp)

        print "Maximal value for s_oc: " + str(max(s_oc.values))
        print "Minimal value for s_oc: " + str(min(s_oc.values))
        print "Maximal value for d_oc: " + str(max(d_oc.values))
        print "Minimal value for d_oc: " + str(min(d_oc.values))
        print "Maximal value for s_txt: " + str(max(s_txt.values))
        print "Minimal value for s_txt: " + str(min(s_txt.values))
        print "Maximal value for d_txt: " + str(max(d_txt.values))
        print "Minimal value for d_txt: " + str(min(d_txt.values))
        print "Maximal value for s_drain: " + str(max(s_drain.values))
        print "Minimal value for s_drain: " + str(min(s_drain.values))
        print "Maximal value for d_drain: " + str(max(d_drain.values))
        print "Minimal value for d_drain: " + str(min(d_drain.values))
        print "Maximal value for temp: " + str(max(temp.values))
        print "Minimal value for temp: " + str(min(temp.values))
        print "Maximal value for tawcd: " + str(max(tawcd.values))
        print "Minimal value for tawcd: " + str(min(tawcd.values))
        print "Maximal value for tawcs: " + str(max(tawcs.values))
        print "Minimal value for tawcs: " + str(min(tawcs.values))        
    
    for icell in range(grassarea.length):
        if (ldebug):
            	fr_txt_cell_grass    = 0.0
	    	fr_oc_cell_grass      = 0.0
                fr_climate_cell_grass = 0.0
                fr_drain_cell_grass   = 0.0
                fr_txt_cell_crop      = 0.0
		fr_oc_cell_crop       = 0.0
                fr_climate_cell_crop  = 0.0
                fr_drain_cell_crop    = 0.0
                fr_txt_cell_nat       = 0.0
		fr_oc_cell_nat        = 0.0
                fr_climate_cell_nat   = 0.0
                fr_drain_cell_nat     = 0.0
        
        Qeff_cell_grass = 0.0
        Qeff_cell_crop  = 0.0
        Qeff_cell_nat   = 0.0
        
        N_denitrification_cell = 0.0
        N_denitrification_cell_arable = 0.0
        N_denitrification_cell_grass = 0.0
        N_denitrification_cell_nat = 0.0           
        total_leaching = 0.0
        Ftmp_cell  = fr_Ndenitrification_temperature(temp.get_data(icell,0.0))
        fQsro_cell = fQsro.get_data(icell,0.0)
        # Only denitrification when pnet is more than params.minimal_soil_infiltration
        pnet_cell  = pnet.get_data(icell,0.0)
        if (pnet_cell < params.minimal_soil_infiltration):
            pnet_cell = 0.0
        tawcd_cell = tawcd.get_data(icell,0.0)
        tawcs_cell = tawcs.get_data(icell,0.0)
        # For grass and crops OC, texture  and drainage
        Foc_cell    = fr_Ndenitrification_oc(s_oc.get_data(icell,0.0))
        Ftxt_cell   = fr_Ndenitrification_text(int(s_txt.get_data(icell,0)))
        Fdrain_cell = fr_Ndenitrification_drainage(int(s_drain.get_data(icell,0)))

        grass_cell = grassarea.get_data(icell,0.0)       
        if (grass_cell > 0.0):
            # There is grass
            # For grass and natural systems only 0.25 of fQsro * pnet is flowing as surface runoff
            Qeff_cell_grass  = (1.0 - params.water_reduction_grass * fQsro_cell) * pnet_cell
            Ftrzd_cell = trzd(Qeff_cell_grass,tawcs_cell)
            fr_Nleach_cell = (1.0 - min(1.0, Ftrzd_cell * Ftmp_cell + Foc_cell + Ftxt_cell + Fdrain_cell))*params.leaching_reduction_grass
            if (ldebug):
                fr_txt_cell_grass     = Ftxt_cell
                fr_oc_cell_grass      = Foc_cell
                fr_climate_cell_grass = Ftrzd_cell * Ftmp_cell
                fr_drain_cell_grass     = Fdrain_cell
            
            N_bal_cell = N_bal_grass.get_data(icell,0.0)
            total_leaching += N_bal_cell * fr_Nleach_cell
            N_denitrification_cell += N_bal_cell * (1.0 - fr_Nleach_cell)
            N_denitrification_cell_grass += N_bal_cell * (1.0 - fr_Nleach_cell)
            N_bal_grass.set_data(icell,N_bal_cell * fr_Nleach_cell)       
            if (N_bal_cell * fr_Nleach_cell < 0.0):
                print "Negative N_bal_grass of ",N_bal_cell * fr_Nleach_cell, N_bal_cell, fr_Nleach_cell

        crop_cell = croparea.get_data(icell,0.0)
        if (crop_cell > 0.0):
            # There is arable land
            Qeff_cell_crop  = (1.0 - params.water_reduction_crop * fQsro_cell) * pnet_cell
            # Because crop must have enough water, we assume that the residence time of the water is maximal 1 year (for example irrigation).
            Ftrzd_cell = min(1.0,trzd(Qeff_cell_crop,tawcs_cell))
            # Added for the oc experiment of crops (30% reduction of oc)
            # Foc_cell    = fr_Ndenitrification_oc(0.7 * s_oc.get_data(icell,0.0))
            fr_Nleach_cell = (1.0 - min(1.0, Ftrzd_cell * Ftmp_cell + Foc_cell + Ftxt_cell + Fdrain_cell))*params.leaching_reduction_crop
            if (ldebug):
                fr_txt_cell_crop     = Ftxt_cell
                fr_oc_cell_crop      = Foc_cell
                fr_climate_cell_crop = Ftrzd_cell * Ftmp_cell            
                fr_drain_cell_crop   = Fdrain_cell
            
            N_bal_cell = N_bal_arable.get_data(icell,0.0)
            total_leaching += N_bal_cell * fr_Nleach_cell
            N_denitrification_cell += N_bal_cell * (1.0 - fr_Nleach_cell)
            N_denitrification_cell_arable += N_bal_cell * (1.0 - fr_Nleach_cell)
            N_bal_arable.set_data(icell,N_bal_cell * fr_Nleach_cell)
            if (N_bal_cell * fr_Nleach_cell < 0.0):
                print "Negative N_bal_arable of ",N_bal_cell * fr_Nleach_cell, N_bal_cell, fr_Nleach_cell        
        nat_cell = natarea.get_data(icell,0.0)
        if (nat_cell > 0.0):
            # There is natural land
            # For grass and natural systems only 0.25 of fQsro * pnet is flowing as surface runoff
            Qeff_cell_nat  = (1.0 - params.water_reduction_natural * fQsro_cell) * pnet_cell
            Ftrzd_cell = trzd(Qeff_cell_nat,tawcd_cell)
            Foc_cell   = fr_Ndenitrification_oc(d_oc.get_data(icell,0.0))
            Ftxt_cell = fr_Ndenitrification_text(int(d_txt.get_data(icell,0)))
            Fdrain_cell = fr_Ndenitrification_drainage(int(d_drain.get_data(icell,0)))
            fr_Nleach_cell = (1.0 - min(1.0, Ftrzd_cell * Ftmp_cell + Foc_cell + Ftxt_cell + Fdrain_cell))*params.leaching_reduction_nat
            if (ldebug):
	        fr_txt_cell_nat     = Ftxt_cell
	        fr_oc_cell_nat      = Foc_cell
                fr_climate_cell_nat = Ftrzd_cell * Ftmp_cell
                fr_drain_cell_nat   = Fdrain_cell
            
            N_bal_cell = N_bal_natural.get_data(icell,0.0)
            total_leaching += N_bal_cell * fr_Nleach_cell
            N_denitrification_cell += N_bal_cell * (1.0 - fr_Nleach_cell)
            N_denitrification_cell_nat += N_bal_cell * (1.0 - fr_Nleach_cell)
            N_bal_natural.set_data(icell,N_bal_cell * fr_Nleach_cell)
            if (N_bal_cell * fr_Nleach_cell < 0.0):
                print "Negative N_bal_natural of ",N_bal_cell * fr_Nleach_cell, N_bal_cell, fr_Nleach_cell
            
        # Put result in raster
        N_denitrification.set_data(icell,N_denitrification_cell)
        N_denitrification_arable.set_data(icell,N_denitrification_cell_arable)
        N_denitrification_grass.set_data(icell,N_denitrification_cell_grass)
        N_denitrification_nat.set_data(icell,N_denitrification_cell_nat)
        Nleach.set_data(icell,total_leaching)

           
        sum_area = nat_cell + crop_cell + grass_cell
        # Calculate the Qeff over all landuse classes
        try:
            Qeff_cell       = (nat_cell/sum_area) * Qeff_cell_nat + \
                              (grass_cell/sum_area) * Qeff_cell_grass + \
                              (crop_cell/sum_area) * Qeff_cell_crop           
        except ZeroDivisionError:
            Qeff_cell = 0.0

        # Put Qeff in a grid 
        Qeff.set_data(icell,Qeff_cell)       

        if (ldebug):   
            try:
                fr_climate_cell = (nat_cell/sum_area) * fr_climate_cell_nat + \
                                  (grass_cell/sum_area) * fr_climate_cell_grass + \
                                  (crop_cell/sum_area) * fr_climate_cell_crop
                fr_oc_cell      = (nat_cell/sum_area) * fr_oc_cell_nat + \
                                  (grass_cell/sum_area) * fr_oc_cell_grass + \
                                  (crop_cell/sum_area) * fr_oc_cell_crop
                fr_txt_cell     = (nat_cell/sum_area) * fr_txt_cell_nat + \
                                  (grass_cell/sum_area) * fr_txt_cell_grass + \
                                  (crop_cell/sum_area) * fr_txt_cell_crop
                fr_drain_cell   = (nat_cell/sum_area) * fr_drain_cell_nat + \
                                  (grass_cell/sum_area) * fr_drain_cell_grass + \
                                  (crop_cell/sum_area) * fr_drain_cell_crop  
            except ZeroDivisionError:
                fr_climate_cell = 0.0
	        fr_oc_cell      = 0.0
	        fr_txt_cell     = 0.0
	        fr_drain_cell   = 0.0

            fr_climate.set_data(icell,fr_climate_cell)
            fr_oc.set_data(icell,fr_oc_cell)
            fr_txt.set_data(icell,fr_txt_cell)
            fr_drain.set_data(icell,fr_drain_cell)
            Ftrzd.set_data(icell,trzd(Qeff_cell,tawcd_cell))
            Ftrzs.set_data(icell,trzd(Qeff_cell,tawcs_cell))
            Ftemp.set_data(icell,Ftmp_cell)
               
        if (N_denitrification_cell < 0.0):
            print "Negative denitrification: ",icell,N_denitrification_cell,Ftrzd_cell,Ftmp_cell,Foc_cell,Ftxt_cell,N_bal_cell
        
        
        
    # Make denitrification as fraction of total load
    fr_denitri  = make_fraction.make_fraction(params,N_denitrification, total_balance)

    # Write fraction to output file  
    fr_denitri.write_ascii_file(os.path.join(params.outputdir,"fr_N_denitr_soil.asc")) 
    
    if (ldebug):
        fr_climate.write_ascii_file(os.path.join(params.outputdir,"fr_N_denitr_soil_climate.asc"))
        fr_oc.write_ascii_file(os.path.join(params.outputdir,"fr_N_denitr_soil_oc.asc"))
        fr_txt.write_ascii_file(os.path.join(params.outputdir,"fr_N_denitr_soil_txt.asc"))
        fr_drain.write_ascii_file(os.path.join(params.outputdir,"fr_N_denitr_soil_drain.asc"))
        Ftrzd.write_ascii_file(os.path.join(params.outputdir,"ftrzd.asc"))
        Ftrzs.write_ascii_file(os.path.join(params.outputdir,"ftrzs.asc"))
        Ftemp.write_ascii_file(os.path.join(params.outputdir,"ftemp.asc"))
    
    # Write Qeff to output file
    Qeff.write_ascii_file(os.path.join(params.outputdir,"Qeff.asc"))
    
    # Write soil denitrification to output file
    N_denitrification.write_ascii_file(os.path.join(params.outputdir,"N_denitr_soil.asc"))
    N_denitrification_arable.write_ascii_file(os.path.join(params.outputdir,"N_denitr_soil_arable.asc"))
    N_denitrification_grass.write_ascii_file(os.path.join(params.outputdir,"N_denitr_soil_grass.asc"))
    N_denitrification_nat.write_ascii_file(os.path.join(params.outputdir,"N_denitr_soil_nat.asc"))
    
    # Write soil leaching to output file
    Nleach.write_ascii_file(os.path.join(params.outputdir,"Nleaching.asc"))
    N_bal_arable.write_ascii_file(os.path.join(params.outputdir,"Nleaching_arable.asc"),output_nodata_value=-1)
    N_bal_grass.write_ascii_file(os.path.join(params.outputdir,"Nleaching_grass.asc"),output_nodata_value=-1)
    N_bal_natural.write_ascii_file(os.path.join(params.outputdir,"Nleaching_nat.asc"),output_nodata_value=-1)
    # Make denitrification as fraction of total load
    fr_Nleach  = make_fraction.make_fraction(params,Nleach, total_balance)
    
    # Delete the data which is not used anymore
    del temp
    del tawcd
    del tawcs
    del s_oc
    del d_oc
    del s_txt

    # 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 soil denitrification." 

    # Aggregation of the rasters
    aggregate.aggregate_grid(basin,N_denitrification,mouth_dict,item="soil_denitrification")
    aggregate.aggregate_grid(basin,N_denitrification_arable,mouth_dict,item="soil_denitrification_arable")
    aggregate.aggregate_grid(basin,N_denitrification_grass,mouth_dict,item="soil_denitrification_grs")
    aggregate.aggregate_grid(basin,N_denitrification_nat,mouth_dict,item="soil_denitrification_nat")
    aggregate.aggregate_grid(basin,fr_denitri,mouth_dict,item="fr_soil_denitrification",method=4,weighing=landarea)
    aggregate.aggregate_grid(basin,Nleach,mouth_dict,item="nleaching")
    aggregate.aggregate_grid(basin,fr_Nleach,mouth_dict,item="fr_nleaching",method=4,weighing=landarea)
    aggregate.aggregate_grid(basin,N_bal_arable,mouth_dict,item="nleaching_arable")
    aggregate.aggregate_grid(basin,N_bal_grass,mouth_dict,item="nleaching_grs")
    aggregate.aggregate_grid(basin,N_bal_natural,mouth_dict,item="nleaching_nat")

    # Calculate the N2O denitrification in the soil
    # Formula is from the IPPC and N2O is a fixed fraction of the leaching.
    # Use the Nleaching grid to calculate the N2O soil denitrification
    N2O = Nleach
    N2O_kgha = ascraster.duplicategrid(landarea)
    for icell in range(N2O.length):
        N2O_cell = Nleach.get_data(icell,0.0) * params.N2O_soil_denitrification_fraction
        N2O.set_data(icell,N2O_cell)
        area_cell = 100.0 * landarea.get_data(icell,0.0) # Conversion from km2 to ha
        try:
            N2O_cell_kgha = N2O_cell/area_cell
        except ZeroDivisionError:
            N2O_cell_kgha = 0.0      
        N2O_kgha.set_data(icell,N2O_cell_kgha)

    N2O.write_ascii_file(os.path.join(params.outputdir,"N2O_denitr_soil.asc"))
    aggregate.aggregate_grid(basin,N2O,mouth_dict,item="N2O_soil_denitrification")
    N2O_kgha.write_ascii_file(os.path.join(params.outputdir,"N2O_denitr_soil_kgha.asc")) 

    # Calculate soil denitrification in kg/ha
    N_kgha = ascraster.duplicategrid(landarea)
    for icell in range(landarea.length):
        area_cell = 100.0 * landarea.get_data(icell,0.0) # Conversion from km2 to ha
        try:
            N_den = N_denitrification.get_data(icell,0.0)/area_cell
        except ZeroDivisionError:
            N_den = 0.0
        N_kgha.set_data(icell,N_den)
    # Write soil denitrification in kg/ha to output file
    N_kgha.write_ascii_file(os.path.join(params.outputdir,"N_denitr_soil_kgha.asc"))

    # Calculate soil denitrification arable in kg/ha
    N_kgha = ascraster.duplicategrid(landarea)
    for icell in range(croparea.length):
        area_cell = 100.0 * croparea.get_data(icell,0.0) # Conversion from km2 to ha
        try:
            N_den = N_denitrification_arable.get_data(icell,0.0)/area_cell
        except ZeroDivisionError:
            N_den = 0.0
        N_kgha.set_data(icell,N_den)
    # Write soil denitrification in kg/ha to output file
    N_kgha.write_ascii_file(os.path.join(params.outputdir,"N_denitr_soil_arable_kgha.asc"))

    # Calculate soil denitrification grass in kg/ha
    N_kgha = ascraster.duplicategrid(landarea)
    for icell in range(grassarea.length):
        area_cell = 100.0 * grassarea.get_data(icell,0.0) # Conversion from km2 to ha
        try:
            N_den = N_denitrification_grass.get_data(icell,0.0)/area_cell
        except ZeroDivisionError:
            N_den = 0.0
        N_kgha.set_data(icell,N_den)
    # Write soil denitrification in kg/ha to output file
    N_kgha.write_ascii_file(os.path.join(params.outputdir,"N_denitr_soil_grass_kgha.asc"))

    # Calculate soil denitrification natural in kg/ha
    N_kgha = ascraster.duplicategrid(landarea)
    for icell in range(natarea.length):
        area_cell = 100.0 * natarea.get_data(icell,0.0) # Conversion from km2 to ha
        try:
            N_den = N_denitrification_nat.get_data(icell,0.0)/area_cell
        except ZeroDivisionError:
            N_den = 0.0
        N_kgha.set_data(icell,N_den)
    # Write soil denitrification in kg/ha to output file
    N_kgha.write_ascii_file(os.path.join(params.outputdir,"N_denitr_soil_nat_kgha.asc"))
 
    return N_denitrification,Qeff
        
def calculate_N_soil_denitrification_riparian(params,mask,\
                       Nsgrw,Qeff,fQgw,mouth_dict,basin):
    '''
    Calculates the N soil dentrification in the riparian zone.
    It returns the ascraster object with the N riparian dentrification.
    The N riparian dentrification is substracted from Nsgrw.
    Here the same method is used as for soil denitrification.
    '''
       
    if params.ldebug: print "Start with calculation of N riparian denitrification."

    # Gerard assumes that there is a DT50 of 2 years, therefore the
    # calibration factor a is equal to log(2)/2 = 0.3466
    a = 0.5 * math.log(2)

    # Read air temperature map
    temp = ascraster.Asciigrid(ascii_file=params.temperature,mask=mask,numtype=float)
    if params.ldebug: print "Air temperature is read in calculation of riparian denitrification."
    
    tawcs = ascraster.Asciigrid(ascii_file=params.tawcs,mask=mask,numtype=float)
    if params.ldebug: print "Tawcs is read in calculation of riparian denitrification."

    s_oc = ascraster.Asciigrid(ascii_file=params.s_oc,mask=mask,numtype=float)
    if params.ldebug: print "S_Oc is read in calculation of riparian denitrification."
    
    s_txt = ascraster.Asciigrid(ascii_file=params.s_txt,mask=mask,numtype=int)
    if params.ldebug: print "S_texture is read in calculation of riparian denitrification."
    
    s_drain = ascraster.Asciigrid(ascii_file=params.s_drain,mask=mask,numtype=int)
    if params.ldebug: print "S_drainage is read in calculation of riparian denitrification."

    s_ph = ascraster.Asciigrid(ascii_file=params.s_ph,mask=mask,numtype=float)
    if params.ldebug: print "S_drainage is read in calculation of fraction riparian denitrification."

    # 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 riparian denitrification." 

    # Grid zero gives the cells with nodata (so no water body). So this map is kept for the whole loop as indicator
    # whether there must be an calculation.
    grid0 = ascraster.Asciigrid(ascii_file=os.path.join(params.glwd_dir,"grid0.asc"),\
                                mask=mask,numtype=float)
 
    # Make a copy (alias not a deepcopy) of the return values
    N_denitrification = temp
    N2O_riparian = s_ph
    N_denitr_kgha = s_oc
    N2O_denitr_kgha = tawcs 
    N_rip_load = ascraster.duplicategrid(temp)
    N_rip_bypass_load = ascraster.duplicategrid(temp)

       
    for icell in range(Nsgrw.length):       
        nsgrw_cell = Nsgrw.get_data(icell,0.0)       
        if (nsgrw_cell > 0.0):
            N_denitrification_cell = 0.0
            Ftmp_cell  = fr_Ndenitrification_temperature(temp.get_data(icell,0.0))
            tawcs_cell = tawcs.get_data(icell,0.0)
            # For grass and crops OC, texture  and drainage
            Foc_cell    = fr_Ndenitrification_oc(s_oc.get_data(icell,0.0))
            Ftxt_cell   = fr_Ndenitrification_text(int(s_txt.get_data(icell,0)))
            Fdrain_cell = fr_Ndenitrification_drainage(int(s_drain.get_data(icell,0)))
            Fph_cell    = fr_Ndenitrification_ph(s_ph.get_data(icell,0.0))

            Qeff_cell  = Qeff.get_data(icell,0.0)
            fQgw_cell  = fQgw.get_data(icell,0.0)
            Ftrzd_cell = params.riparian_layer_depth * trzd(Qeff_cell*(1.0-fQgw_cell),tawcs_cell)
            fr_denitr_cell = min(1.0, Ftrzd_cell * Ftmp_cell + Foc_cell + Ftxt_cell + Fdrain_cell)

            # Only the fraction of the cell with no water body has riparian denitrification
            fr_non_water_area = grid0.get_data(icell,1.0)
            # fr_denitr_cell =  fr_denitr_cell * fr_non_water_area

            # pH effect on denitrification
            ph = s_ph.get_data(icell,0.0)
            if (ph > 0.0): 
                fr_denitr_cell *=  interpolate(ph,ph1_xlist,ph1_ylist)

            # Put result in raster
            N_riparian_cell = nsgrw_cell * fr_non_water_area
            N_rip_bypass_load_cell = nsgrw_cell * (1.0 - fr_non_water_area)
            N_rip_load.set_data(icell,N_riparian_cell)
            N_rip_bypass_load.set_data(icell,N_rip_bypass_load_cell)
            N_rip_denitr_cell = N_riparian_cell * fr_denitr_cell
            N2O_riparian_cell = N_rip_denitr_cell *  Fph_cell * (1.0 - fr_denitr_cell) 
            # print "N2O ",icell,nsgrw_cell * fr_denitr_cell,(1.0 - fr_denitr_cell),N2O_riparian_cell
            N2O_riparian.set_data(icell,N2O_riparian_cell)
            N_denitrification.set_data(icell,N_rip_denitr_cell)
            Nsgrw.set_data(icell,nsgrw_cell - N_rip_denitr_cell)

            # Calculate the load in kg/ha
            area_cell = 100.0 * landarea.get_data(icell,0.0) # Conversion from km2 to ha
            try:
                N2O_cell_kgha = N2O_riparian_cell/area_cell
            except ZeroDivisionError:
                N2O_cell_kgha = 0.0      
            N2O_denitr_kgha.set_data(icell,N2O_cell_kgha)
            try:
                N_cell_kgha = N_riparian_cell/area_cell
            except ZeroDivisionError:
                N_cell_kgha = 0.0      
            N_denitr_kgha.set_data(icell,N_cell_kgha)    

        else:
            # Put result in raster    
            N_denitrification.set_data(icell,0.0)
            N2O_riparian.set_data(icell,0.0)
            N_denitr_kgha.set_data(icell,0.0)
            N2O_denitr_kgha.set_data(icell,0.0)
            N_rip_load.set_data(icell,0.0)
            N_rip_bypass_load.set_data(icell,0.0)      
   
    # Write soil denitrification to output file
    N_denitrification.write_ascii_file(os.path.join(params.outputdir,"N_denitr_riparian.asc"))
    N2O_riparian.write_ascii_file(os.path.join(params.outputdir,"N2O_denitr_riparian.asc"))
    N_denitr_kgha.write_ascii_file(os.path.join(params.outputdir,"N_denitr_riparian_kgha.asc"))
    N2O_denitr_kgha.write_ascii_file(os.path.join(params.outputdir,"N2O_denitr_riparian_kgha.asc"))
    
    # Write Nsgrw without riparian denitrification to output file
    Nsgrw.write_ascii_file(os.path.join(params.outputdir,"Nsgrw.asc"))
    
    # Delete the data which is not used anymore
    del temp
    del tawcs
    del s_oc
    del s_txt

    # Aggregation of the rasters
    aggregate.aggregate_grid(basin,N_rip_load,mouth_dict,item="n_rip_load")
    aggregate.aggregate_grid(basin,N_rip_bypass_load,mouth_dict,item="n_rip_bypass_load")
    aggregate.aggregate_grid(basin,N_denitrification,mouth_dict,item="n_denitr_rip")
    aggregate.aggregate_grid(basin,N2O_riparian,mouth_dict,item="n2o_denitr_rip")
    aggregate.aggregate_grid(basin,Nsgrw,mouth_dict,item="nsgrw_end")
 
    return N_denitrification,Nsgrw    
