# ******************************************************
## Revision "$LastChangedDate: 2012-05-24 18:36:08 +0200 (Thu, 24 May 2012) $"
## Date "$LastChangedRevision: 222 $"
## Author "$LastChangedBy: arthurbeusen $"
## URL "$HeadURL: https://pbl.sliksvn.com/globalnutrients/P_model/trunk/tools/pbackground.py $"
# ******************************************************

import os
import math
import ascraster
import my_sys
import aggregate
import litho_class
from error import *

# The average Temp of Japan is 11 degree Celsius.
Temp_Japan = 284.0

def fr_temperature(temp,Ea):
    '''
    For the 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  = constant = 1
    Ea = activation energy 
    R  = universal gasconstant 8.3144 J/(K.MOL)
    K  = temperature in Kelvin (Celcius+273)
    As for the Activiation energy: the global chemical weathering rates are relative to Japan (11 degree Celsius).
    therefor you need to calculate: exp(-Ea/R  *(1/Tlocal - 1/TJapan)) or [exp (-Ea *1/R * 1/Tlocal) / exp (-Ea * 1/R * 1/Tjapan) ]
    The average Temp of Japan is 11 degree Celsius.
    '''
    return math.exp((-Ea/8.3144)*(1.0/(temp+273) - 1.0/(Temp_Japan)))   

def calculate_P(params,mask,pnet,mouth_dict,basin,landarea):
    '''
    Calculate the background load of P groundwater to surface water.
    '''
       
    if params.ldebug: print "Start with calculation of P weathering."

    # Make the soil shielding factor for each soil type
    shielding_factor = [1.0, 0.1, 0.1, 0.1, 0.1] 

    # Make default grid
    shielding_grid = ascraster.duplicategrid(landarea) 
    shielding_grid.add_values(shielding_grid.length*[0.0])

    # Loop over all the classes of soil shielding:
    for iclass in range(0,5):
        # Check header of input file
        equal = ascraster.compare_grids(params.landarea,os.path.join(params.soil_shield_dir,"grid" + str(iclass) + ".asc"))
        if (not equal):
            raise MyError("Soil Shielding maps don't have same header as other input files.")


        grid = ascraster.Asciigrid(ascii_file=os.path.join(params.soil_shield_dir,"grid" + str(iclass) + ".asc"),\
                                   mask=mask,numtype=float)
        for icell in xrange(grid.length):
            grid_cell = grid.get_data(icell,0.0)
            shield_cell = shielding_grid.get_data(icell)
            shielding_grid.set_data(icell,shield_cell + shielding_factor[iclass]*grid_cell)
                

    temp_grid = ascraster.Asciigrid(ascii_file=params.temperature,mask=mask,numtype=float)

    # Read litho map and determine the background concentration due to weathering
    litho_values = litho_class.Litho().make_value_list(params,mask,["backgroundPO4","Ea"])
          
    backgroundP = ascraster.duplicategrid(landarea)
    for icell in range(backgroundP.length):
        pnet_cell = pnet.get_data(icell,0.0)
        if (pnet_cell > 0.0):
	    # BackgroundPO4 is given in mgP/l
            concP = litho_values[icell][0]
            area = landarea.get_data(icell,0.0)
            # pnet is in mm/yr  and BackgroundPO4 is given in mgP/l. We want a load in kg P.
            # pnet * concP = mm/yr * mg/l = 1e-3 m/yr * 1e-6 kg/1e-3 m3 = 1e-9 kg/1e-3 m2 = 1e-9 kg/1e-9 km2 = kg/km2 
            backgnd = concP * pnet_cell * area
            # If there is no temperature, use the Japan average temperature, so correction factor will be 1.0
            temp = temp_grid.get_data(icell,Temp_Japan-273.0)
            temp_correction = fr_temperature(temp,litho_values[icell][1])
            shield_correction = shielding_grid.get_data(icell)
            backgroundP.set_data(icell,backgnd*temp_correction*shield_correction)
        else:
            backgroundP.set_data(icell,0.0)


    # Aggregate rasters
    aggregate.aggregate_grid(basin,backgroundP,mouth_dict,item="pbackground")

    # Write to output file
    backgroundP.write_ascii_file(os.path.join(params.outputdir,"pbackground.asc"))
 
    return backgroundP
