#"""
#   Simulations on the GBR with multirate explicit Runge-Kutta time-stepping.
#   All the parameters can be found and modified in the file "InpuGBR.py"
#"""

#   A) HOW TO LAUNCH A SIMULATION
#       Simulation on 1 processor  : rundgpy MultirateGBR.py
#       Simulation on n processors : mpirun -np n rundgpy MultirateGBR.py   ( !!! You need MPI !!! )
#
#
#   B) TO START SIMULATION FROM A PRE-EXISTING SOLUTION:
#     1. Set the following in Input.py:
#         -> set flag startFromExistingSolution = 1
#         -> set solutionFileName
#         -> set stepNbOfInitialSolution
#         -> set old_dt
#         -> then, Ti = old Ti + stepNbOfInitialSolution*old_dt
#     2. Set the following in MultirateGBR.py:
#         -> Nothing
#
#         Notes:  Must point to a 'merged' .idx file (use generateIdxFiles.py, and make sure you adjust nb_parts)
#                 Make sure idx file is complete (ie open it and check it looks right) - especially if previous job crashed out!
#"""

# Import for python
from dgpy import *
from gmshpy import *
from math import *
from ProjectMesh import *
from InputGBR_oldProj import *
from termcolor import colored
import time, os, sys

if(Msg.GetCommRank() == 0):
  print''
  print(colored('SLIM: Multirate GBR simulation starting. Loading 1:Mesh, 2:Bathymetry 3:Reef Map 4:Conservation law.',"yellow"))
  print''

# Lauch Makefile to generate meshes and download forcings
if(Msg.GetCommRank() == 0):
  try : os.mkdir('./Meshes');
  except: 0;
  try : os.mkdir('./Data');
  except: 0;
  try : os.mkdir('./GBR_Functions');
  except: 0;
  if(not os.path.exists("./Meshes/"+filename+".msh")):
    try : os.system("make "+filename+".msh");
    except: 0;
  if (gbrWind == 2):
    try : os.system("make forcings_windNCEP")
    except: 0;
  elif (gbrWind == 3):
    try : os.system("make forcings_windCSFR")
    except: 0;
  try : os.system("make forcings");
  except: 0;
  try : os.mkdir(outputDir);
  except: 0;
  try : os.mkdir(outputDir+'/Evaluator');
  except: 0;

Msg.Barrier()

# Definition of Python functions
def printTime (t) :
  return "%dh%02d'%02d\"" % (t/3600,  (t%3600)/60,  t%60)

def printDate (t):
  date=time.gmtime(t)
  return "%02d/%02d/%04d %02dh%02d'%02d\""%(date[2], date[1], date[0], date[3], date[4], date[5])

# Compile libraries
if (Msg.GetCommRank() == 0):
  functionC.buildLibraryFromFile ("Libraries/GBR.cc", "Libraries/lib_gbr.so")
glib = "Libraries/lib_gbr.so"
Msg.Barrier()

# Load model
model = GModel()
if(Msg.GetCommSize() == 1):
  print(colored('Loading mesh for 1 processor (unpartitioned) ...',"red"))
  if(not os.path.exists('./Meshes/'+filename+'_tan.msh')):
    Msg.Info("Projected mesh doesn't exist. Projecting it.")
    m = GModel()
    m.load("./Meshes/"+filename+".msh")
    projectMeshPlane(m)
    m.save("./Meshes/"+filename+"_tan"+".msh")
  model.load('./Meshes/'+filename+'_tan.msh')  
else:
  if(os.path.exists("./Meshes/"+filename+"_tan_partitioned_"+str(Msg.GetCommSize())+".msh")):
    if(Msg.GetCommRank() == 0):
      print'Loading partitioned mesh: ./Meshes/'+filename+'_partitioned_'+str(Msg.GetCommSize())+'.msh'
    model.load('./Meshes/'+filename+'_tan_partitioned_'+str(Msg.GetCommSize())+'.msh')
  else:
    Msg.Fatal('No valid partitioned mesh file: ./Meshes/'+filename+'_partitioned_[nbProcs].msh ; run PartitionGBR_Standalone.py first.')

# Coordinate system
XYZ = function.getCoordinates()
lonLatDegrees = functionC(glib,"lonLatDegrees",3,[XYZ])
lonLat = functionC(glib,"lonLat",3,[XYZ])

latLonDegrees = functionC(glib,"latLonDegrees",3,[XYZ])
latLon = functionC(glib,"latLon",3,[XYZ])

groups = dgGroupCollection(model, dimension, order)
#groups.buildGroupsOfInterfaces()
bathDC = dgDofContainer(groups, 1)


# ****************** Bathymetry ******************

if(os.path.exists("./Bath/"+filename+"_bath_smooth/"+filename+"_bath_smooth.idx")):
  Msg.Info("Smoothed bathymetry already exists.")
  bathDC.importIdx("./Bath/"+filename+"_bath_smooth/"+filename+"_bath_smooth.idx")
  Msg.Info("Smoothed bathymetry read.")
else:
  Msg.Fatal("You need to run DiffuseBathGBR_Standalone first!")
#diffuseBathymetryStereo(bathDC, groups, stereoToLonLatDegrees)
#if using full.bin: diffuseBathymetry(bathDC, groups, lonLat)

# ************************************************

# ****************** Reefs ******************
if (Msg.GetCommRank() == 0):
  print(colored('Loading reef map ...',"red"))

reefsDC = dgDofContainer(groups,1)

reefsDCfilename = "GBR_Functions/reefsFunction_"+filename+"/reefsFunction_"+filename+".idx"
if(os.path.exists(reefsDCfilename)):
  if (Msg.GetCommRank() == 0):
    print 'Reef Map file found: '+reefsDCfilename+' ; loading this.'
  reefsDC.importIdx(reefsDCfilename)
else:
  if (Msg.GetCommRank() == 0):
    print 'Reef Map file '+reefsDCfilename+' not found. Loading reef map ...'
  #1. Create a reefs object (initialise object and create reefs function)
  reefsObj = reefsSmooth(reefMapFile, lonLatDegrees) 
  #2. Interpolate the function onto DC
  reefsDC.interpolate(reefsObj)
  exporterReefsDC=dgIdxExporter(reefsDC, 'GBR_Functions/reefsFunction_'+filename)
  exporterReefsDC.exportIdx(0,0)
  #3. Delete reefs object data
  reefsObj.clearReefsData()

# ************************************************

t = Ti
timeFunction = functionConstant(t)

# ****************** Forcings ******************
if (Msg.GetCommRank() == 0):
  print(colored('Setting up forcings ...',"red"))
# .................. Tides    ..................

# Tpxo tide 
tideEta = slimFunctionTpxo("./Data/h_tpxo7.2.nc","ha","hp","lon_z","lat_z")
tideEta.setCoordAndTimeFunctions(lonLatDegrees,timeFunction)
tideU = slimFunctionTpxo("./Data/u_tpxo7.2.nc","Ua","up","lon_u","lat_u")
tideU.setCoordAndTimeFunctions(lonLatDegrees,timeFunction)
tideV = slimFunctionTpxo("./Data/u_tpxo7.2.nc","Va","vp","lon_v","lat_v")
tideV.setCoordAndTimeFunctions(lonLatDegrees,timeFunction)
tideUV = functionC(glib,"latLonVector",3,[latLon,tideV,tideU])

# .................. Wind     ..................

windInterpolatorU = slimStructDataInterpolatorMultilinear()
windInterpolatorV = slimStructDataInterpolatorMultilinear()
windUcontainer = None
windVcontainer = None
if gbrWind == 1:
  windUcontainer = slimStructDataContainerTemporalSerie("./Data/AIMS_wind_Davies_U.txt")
  windVcontainer = slimStructDataContainerTemporalSerie("./Data/AIMS_wind_Davies_V.txt")
elif gbrWind == 2:
  windUcontainer = slimStructDataContainerNetcdf("./Data/uwind-2000-2011.nc", "uwnd", "lat", "lon")
  windVcontainer = slimStructDataContainerNetcdf("./Data/vwind-2000-2011.nc", "vwnd", "lat", "lon")
elif gbrWind == 3:
  windUcontainer = slimStructDataContainerNetcdf("./Data/wnd10m_csfr_2008-2010.grb2.nc", "U_GRD_L103", "lat", "lon")
  windVcontainer = slimStructDataContainerNetcdf("./Data/wnd10m_csfr_2008-2010.grb2.nc", "V_GRD_L103", "lat", "lon")
windU = slimFunctionStructData(windUcontainer, windInterpolatorU, latLonDegrees)
windV = slimFunctionStructData(windVcontainer, windInterpolatorV, latLonDegrees)

windU.setTimeFunction(timeFunction)
windV.setTimeFunction(timeFunction)

#windU_PostSetTimeFn_DC = dgDofContainer(groups, 1)
#windU_PostSetTimeFn_DC.interpolate(windU)
#windV_PostSetTimeFn_DC = dgDofContainer(groups, 1)
#windV_PostSetTimeFn_DC.interpolate(windV)
#windU_PostSetTimeFn_DC.exportMsh('windU_PostSetTimeFn')
#windV_PostSetTimeFn_DC.exportMsh('windV_PostSetTimeFn')

windUV = functionC(glib,"latLonVector",3,[latLon,windV,windU])

# ................ Bluelink     .................
blinkUV = None
blinkEta_cutOff = None
if gbrResidual == 2:
  blinkInterpolator = slimStructDataInterpolatorMultilinear()
  blinkUContainer = slimStructDataContainerNetcdf("./Data/bluelink_transport_27nov2007_35days.nc", "u", "lat", "lon")
  blinkVContainer = slimStructDataContainerNetcdf("./Data/bluelink_transport_27nov2007_35days.nc", "v", "lat", "lon")
  blinkEtaContainer = slimStructDataContainerNetcdf("./Data/bluelink_eta_27nov2007_35days.nc", "eta", "lat", "lon") 
  blinkU = slimFunctionStructData(blinkUContainer, blinkInterpolator, latLonDegrees)
  blinkV = slimFunctionStructData(blinkVContainer, blinkInterpolator, latLonDegrees)
  blinkEta = slimFunctionStructData(blinkEtaContainer, blinkInterpolator, latLonDegrees)
  blinkU.setTimeFunction(timeFunction)
  blinkV.setTimeFunction(timeFunction)
  blinkEta.setTimeFunction(timeFunction)
  #blinkEta.setMinValue(-3.0)
  #blinkEta.setMaxValue(3.0)
  blinkUV = functionC(glib, "latLonVector",3,[latLon,blinkV,blinkU])

  minFC = functionConstant(-3.0)
  blinkEta_cutOff = functionC(glib, "cutOffMin", 1, [blinkEta, minFC])

# .......... GBROOS Boundary Currents    ..........
if gbrResidual == 3:
  boundaryInterpolator = slimStructDataInterpolatorMultilinear()
  boundaryElevContainer_MYR = slimStructDataContainerTemporalSerie("./Data/currentsForcing_MYR_CCH/MYR_elev_novdec07janfeb08.txt")
  boundaryUcontainer_MYR = slimStructDataContainerTemporalSerie("./Data/currentsForcing_MYR_CCH/MYR_ucur_novdec07janfeb08.txt")
  boundaryVcontainer_MYR = slimStructDataContainerTemporalSerie("./Data/currentsForcing_MYR_CCH/MYR_vcur_novdec07janfeb08.txt")
  boundaryElev_MYR = slimFunctionStructData(boundaryElevContainer_MYR, boundaryInterpolator, latLonDegrees)
  boundaryU_MYR = slimFunctionStructData(boundaryUcontainer_MYR, boundaryInterpolator, latLonDegrees)
  boundaryV_MYR = slimFunctionStructData(boundaryVcontainer_MYR, boundaryInterpolator, latLonDegrees)
  boundaryElev_MYR.setTimeFunction(timeFunctionGBROOS)
  boundaryU_MYR.setTimeFunction(timeFunctionGBROOS)
  boundaryV_MYR.setTimeFunction(timeFunctionGBROOS)
  boundaryUV_MYR = functionC(glib, "latLonVector", 3, [latLon, boundaryV_MYR, boundaryU_MYR])
  #boundaryUV_modulated_MYR = functionC(glib,"currentsAlongVecDotNormals", 3, [function.getNormals(), data_Central, boundaryUV_MYR])
  boundaryUVElev_MYR = functionC(glib, "merge", 3, [boundaryElev_MYR, boundaryUV_MYR])

  boundaryElevContainer_CCH = slimStructDataContainerTemporalSerie("./Data/currentsForcing_MYR_CCH/CCH_elev_novdec07janfeb08.txt")
  boundaryUcontainer_CCH = slimStructDataContainerTemporalSerie("./Data/currentsForcing_MYR_CCH/CCH_ucur_novdec07janfeb08.txt")
  boundaryVcontainer_CCH = slimStructDataContainerTemporalSerie("./Data/currentsForcing_MYR_CCH/CCH_vcur_novdec07janfeb08.txt")
  boundaryElev_CCH = slimFunctionStructData(boundaryElevContainer_CCH, boundaryInterpolator, latLonDegrees)
  boundaryU_CCH = slimFunctionStructData(boundaryUcontainer_CCH, boundaryInterpolator, latLonDegrees)
  boundaryV_CCH = slimFunctionStructData(boundaryVcontainer_CCH, boundaryInterpolator, latLonDegrees)
  boundaryElev_CCH.setTimeFunction(timeFunctionGBROOS)
  boundaryU_CCH.setTimeFunction(timeFunctionGBROOS)
  boundaryV_CCH.setTimeFunction(timeFunctionGBROOS)
  boundaryUV_CCH = functionC(glib,"latLonVector", 3, [latLon, boundaryV_CCH, boundaryU_CCH])
  #boundaryUV_modulated_CCH = functionC(glib,"currentsAlongVecDotNormals", 3, [function.getNormals(), data_South, boundaryUV_CCH])
  boundaryUVElev_CCH = functionC(glib, "merge", 3, [boundaryElev_CCH, boundaryUV_CCH])

# ............  Residual currents    ............

#Put data into functions
data_A = functionConstant([F_A, vx_A, vy_A])
data_B = functionConstant([F_B, vx_B, vy_B])
data_C = functionConstant([F_C, vx_C, vy_C])

# NCJ Current forcing
bathTimesVector_A = functionC(glib, "bathVec", 3, [bathDC.getFunction(), function.getNormals(), data_A])
bathTimesVector_B = functionC(glib, "bathVec", 3, [bathDC.getFunction(), function.getNormals(), data_B])
bathTimesVector_C = functionC(glib, "bathVec", 3, [bathDC.getFunction(), function.getNormals(), data_C])

bathIntegralNorthFM = fullMatrixDouble(1,1)
bathIntegralCentralFM = fullMatrixDouble(1,1)
bathIntegralSouthFM = fullMatrixDouble(1,1)

dgFunctionIntegratorInterface(groups, bathTimesVector_A).compute('Open Sea North', bathIntegralNorthFM)
dgFunctionIntegratorInterface(groups, bathTimesVector_B).compute('Open Sea Central', bathIntegralCentralFM)
dgFunctionIntegratorInterface(groups, bathTimesVector_C).compute('Open Sea South', bathIntegralSouthFM)
bathIntegralsFuncNorth = functionConstant(bathIntegralNorthFM(0,0))
bathIntegralsFuncCentral = functionConstant(bathIntegralCentralFM(0,0))
bathIntegralsFuncSouth = functionConstant(bathIntegralSouthFM(0,0))

# Generate SEC UV functions
currentUV_North = functionC(glib, "SEC_UV", 3, [bathIntegralsFuncNorth, function.getNormals(), data_A])
currentUV_Central = functionC(glib, "SEC_UV", 3, [bathIntegralsFuncCentral, function.getNormals(), data_B])
currentUV_South = functionC(glib, "SEC_UV", 3, [bathIntegralsFuncSouth, function.getNormals(), data_C])

if (Msg.GetCommRank() == 0):
  outputVec = open('GBR_Functions/CurrentUnitVectors.pos', 'w')
  outputVec.write('View \"CurrentUnitVectors\" {\n')
  outputVec.write('VP (' + '%.5f'%0 + ',' + '%.5f'%0 + ',0) {' + '%.5f'%vx_A + ',' + '%.5f'%vy_A + ',0};\n')
  outputVec.write('VP (' + '%.5f'%0 + ',' + '%.5f'%0 + ',0) {' + '%.5f'%vx_B + ',' + '%.5f'%vy_B + ',0};\n')
  outputVec.write('VP (' + '%.5f'%0 + ',' + '%.5f'%0 + ',0) {' + '%.5f'%vx_C + ',' + '%.5f'%vy_C + ',0};\n')
  outputVec.write('};')
  outputVec.close()

# **********************************************


# ****************** Cons Law ******************
if (Msg.GetCommRank() == 0):
  print(colored('Setting up conservation law ...',"red"))

claw = dgConservationLawShallowWater2d()
solution = dgDofContainer(groups, claw.getNbFields())

# **********************************************
# ** OPTIONAL 1/3: Load an existing DC as the initial solution DC

if (startFromExistingSolution == 1):
  solution.importIdx( solutionFileName )

# **********************************************

ld = functionConstant(0)

#Cd = functionC(glib, "bottomDrag_old", 1, [solution.getFunction(), bathDC.getFunction()])
Cd = functionC(glib, "bottomDrag", 1, [ solution.getFunction(), bathDC.getFunction(), reefsDC.getFunction() ]) #*** 1/5
if(gbrDiff == 0):
  Di = functionConstant(diff)
elif(gbrDiff == 1):
  Di = functionC(glib, "smagorinsky", 1, [solution.getFunctionGradient()]) #***
Coriolis = functionC(glib, "coriolisLonLatDegrees", 1, [lonLatDegrees])
if(gbrWind == 0):
  source = functionConstant([0.0, 0.0])
else:
  source = functionC(glib, "windStressSmithBanke", 2, [windUV, solution.getFunction(), bathDC.getFunction()]) #***

claw.setCoriolisFactor(Coriolis)
claw.setQuadraticDissipation(Cd)
claw.setLinearDissipation(ld)
claw.setDiffusivity(Di)
claw.setSource(source)
claw.setBathymetry(bathDC.getFunction())
claw.setBathymetryGradient(bathDC.getFunctionGradient())

tideuv = functionC(glib,"transport2velocity",3,[tideUV, bathDC.getFunction(), solution.getFunction()]) #***
tpxoTide = functionC(glib, "merge", 3, [tideEta, tideuv])

blinkAndTides_VelAndElev = None
if gbrResidual == 2:
  blinkuv = functionC(glib,"transport2velocity",3,[blinkUV, bathDC.getFunction(), solution.getFunction()])
  blink_VelAndElev = functionC(glib, "merge", 3, [blinkEta_cutOff, blinkuv])
  blinkAndTides_VelAndElev = functionC(glib,"merge_3DFunctions",3,[tpxoTide, blink_VelAndElev])

CurrentsAndTideNorth = functionC(glib,"merge3",3,[tideEta, tideuv, currentUV_North])
CurrentsAndTideCentral = functionC(glib,"merge3",3,[tideEta, tideuv, currentUV_Central])
CurrentsAndTideSouth = functionC(glib,"merge3",3,[tideEta, tideuv, currentUV_South])

if(gbrResidual == 0):
  if(gbrTide == 0):
    claw.addBoundaryCondition('Open Sea North',claw.newBoundaryWall())
    claw.addBoundaryCondition('Open Sea Central',claw.newBoundaryWall())
    claw.addBoundaryCondition('Open Sea South',claw.newBoundaryWall())
  elif(gbrTide == 1):
    claw.addBoundaryCondition('Open Sea North',claw.newOutsideValueBoundary("Surface", tpxoTide))
    claw.addBoundaryCondition('Open Sea Central',claw.newOutsideValueBoundary("Surface", tpxoTide))
    claw.addBoundaryCondition('Open Sea South',claw.newOutsideValueBoundary("Surface", tpxoTide))
elif(gbrResidual == 1):
  if(gbrTide == 0):
    claw.addBoundaryCondition('Open Sea North',claw.newOutsideValueBoundary("Surface", currentUV_North))
    claw.addBoundaryCondition('Open Sea Central',claw.newOutsideValueBoundary("Surface", currentUV_Central))
    claw.addBoundaryCondition('Open Sea South',claw.newOutsideValueBoundary("Surface", currentUV_South))
  elif(gbrTide == 1):
    claw.addBoundaryCondition('Open Sea North', claw.newOutsideValueBoundary("Surface", CurrentsAndTideNorth))
    claw.addBoundaryCondition('Open Sea Central', claw.newOutsideValueBoundary("Surface", CurrentsAndTideCentral))
    claw.addBoundaryCondition('Open Sea South', claw.newOutsideValueBoundary("Surface", CurrentsAndTideSouth))
elif(gbrResidual == 2):
  if(gbrTide == 0):
    claw.addBoundaryCondition('Open Sea North',claw.newOutsideValueBoundary("Surface", blink_VelAndElev))
    claw.addBoundaryCondition('Open Sea Central',claw.newOutsideValueBoundary("Surface", blink_VelAndElev))
    claw.addBoundaryCondition('Open Sea South',claw.newOutsideValueBoundary("Surface", blink_VelAndElev))
  elif(gbrTide == 1):
    claw.addBoundaryCondition('Open Sea North', claw.newOutsideValueBoundary("Surface", blinkAndTides_VelAndElev))
    claw.addBoundaryCondition('Open Sea Central', claw.newOutsideValueBoundary("Surface", blinkAndTides_VelAndElev))
    claw.addBoundaryCondition('Open Sea South', claw.newOutsideValueBoundary("Surface", blinkAndTides_VelAndElev))
elif(gbrResidual == 3):
  if(gbrTide == 0):
    claw.addBoundaryCondition('Open Sea North', claw.newBoundaryWall())
    claw.addBoundaryCondition('Open Sea Central', claw.newOutsideValueBoundary("Surface", boundaryUVElev_MYR))
    claw.addBoundaryCondition('Open Sea Neutral', claw.newBoundaryWall())
    claw.addBoundaryCondition('Open Sea South', claw.newOutsideValueBoundary("Surface", boundaryUVElev_CCH))
  elif(gbrTide == 1):
    claw.addBoundaryCondition('Open Sea North', claw.newOutsideValueBoundary("Surface", tpxoTide))
    claw.addBoundaryCondition('Open Sea Central', claw.newOutsideValueBoundary("Surface", boundaryUVElev_MYR))
    claw.addBoundaryCondition('Open Sea Neutral', claw.newOutsideValueBoundary("Surface", tpxoTide))
    claw.addBoundaryCondition('Open Sea South', claw.newOutsideValueBoundary("Surface", boundaryUVElev_CCH))


claw.addBoundaryCondition('Coast',claw.newBoundaryWall())
claw.addBoundaryCondition('Island',claw.newBoundaryWall())
claw.addBoundaryCondition('Islands',claw.newBoundaryWall())

# **********************************************
#NB: Must make 1st argument of Cd solution.getFunction() to use the following:
bottomDragDC = dgDofContainer(groups, 1)
bottomDragDC.interpolate(Cd)
exporterQuadDissip = dgIdxExporter(bottomDragDC, 'GBR_Functions/quadDissipation_'+filename)
exporterQuadDissip.exportIdx(0,0)

# ****************** Time Int ******************
#WARNING: in MultiRate, all DofContainers are split into groups
#       ---> Must add every DC used by claw to splitForMultiRate!!!
if (Msg.GetCommRank() == 0):
  print(colored('Setting up time integrator ...',"red"))

rk = dgMultirateERK(groups, claw, RK_TYPE)
dt = rk.splitGroupsForMultirate(mL, solution, [solution, bathDC, reefsDC])
rk.printMultirateInfo(0)

#OPTIONAL 2/3: if starting from existing solution, set dt = old_dt
if (startFromExistingSolution == 1):
  dt = old_dt
# **********************************************

##Write dissipation term
#dissipationTermFC = functionC(glib, "dissipationTerm", 2, [solution.getFunction(), solution.getFunctionGradient(), Di, bathDC.getFunction(), bathDC.getFunctionGradient()])
#dissipationTermDC = dgDofContainer(groups, 2)
#dissipationTermDC.interpolate(dissipationTermFC)
#exporterDissipationTerm = dgIdxExporter(dissipationTermDC, outputDir+'/DissipationTerm/dissipationTerm')
#exporterDissipationTerm.exportIdx(0,0)

#Write windUV term
windUVDC = None
if (gbrWind != 0):
  windUVDC = dgDofContainer(groups, 3)
  windUVDC.interpolate(windUV)
  exporterWindUV = dgIdxExporter(windUVDC, outputDir+'/windUV')
  exporterWindUV.exportIdx(0, t)
  
#Write windStress term
windSourceDC = None
if (gbrWind != 0):
  windSourceDC = dgDofContainer(groups, 2)
  windSourceDC.interpolate(source)
  exporterWindSourceAcc = dgIdxExporter(windSourceDC, outputDir+'/windSourceAcc')
  exporterWindSourceAcc.exportIdx(0, t)
  
#Write blink term
blinkSolDC = None
if(gbrResidual == 2):
  blinkSolDC = dgDofContainer(groups, 3)
  blinkSolDC.interpolate(blink_VelAndElev)
  exporterBlinkSol=dgIdxExporter(blinkSolDC, outputDir+'/blinkSol')
  exporterBlinkSol.exportIdx(0, t)
  
#Write GBROOS boundary currents/elev term
boundaryGBROOS_MYRSol_DC = None
boundaryGBROOS_CCHSol_DC = None
if(gbrResidual == 3):
  boundaryGBROOS_MYRSol_DC = dgDofContainer(groups, 3)
  boundaryGBROOS_MYRSol_DC.interpolate(boundaryUVElev_MYR)
  exporterGBROOS_MYRSol=dgIdxExporter(boundaryGBROOS_MYRSol_DC, outputDir+'/boundaryGBROOS_MYRSol')
  exporterGBROOS_MYRSol.exportIdx(0, t)
  boundaryGBROOS_CCHSol_DC = dgDofContainer(groups, 3)
  boundaryGBROOS_CCHSol_DC.interpolate(boundaryUVElev_CCH)
  exporterGBROOS_CCHSol=dgIdxExporter(boundaryGBROOS_CCHSol_DC, outputDir+'/boundaryGBROOS_CCHSol')
  exporterGBROOS_CCHSol.exportIdx(0, t)

currentFunction = functionC(glib,"current",3,[solution.getFunction()]) #***
tideCurrentFunction = functionC(glib,"current",3,[tpxoTide])

nbSteps = int(ceil((Tf-Ti)/dt))
nbExport = 0

if(Msg.GetCommRank() == 0):
  if(os.path.exists(outputDir+'/mr-gbr')):
    try : os.system("rm -r "+outputDir+'/mr-gbr');
    except: 0;
exporterSol=dgIdxExporter(solution, outputDir+'/mr-gbr')
#exporterSol=dgIdxExporter(currentFunction, outputDir+'/mr-gbr-func')

norm=solution.norm()

if (Msg.GetCommRank() == 0):
  print ''
  print(colored('******************** Initialising ********************', "blue"))
  print ''
  print '     - Number of processors:', Msg.GetCommSize()
  print '     - Multirate Runge Kutta Scheme:', RK_TYPE
  print '     - Initial Time (Ti):',  printDate(Ti)
  print '     - Final Time (Tf):',  printDate(Tf)
  print '     - Simulation Length (days):',  (Tf-Ti)/3600.0/24.0  
  print '     - Reference Time Step:', printTime(dt), ', which is', str(dt), 's'
  print '     - Number of Iterations:', nbSteps
  print '     - Theoretical Speedup:', rk.speedUp()
  print '     - Norm of Solution at Ti:',  norm
  print ''
  print '     - Tides:', gbrTideLabels[gbrTide], ' || Wind:', gbrWindLabels[gbrWind], ' || Residual current:', gbrResidualLabels[gbrResidual]
  print ''
  print '     - Output Directory:', outputDir
  print ''
  print(colored('******************************************************', "blue"))
  print ''
  
  #Write info file for LPTracker
  simInfo = [str(dt)+'\n', str(nbSteps)+'\n', str(export)+'\n', str(Ti)+'\n', str(outputDir)+'\n', str(filename)]
  simInfoFile = open(outputDir+'/simulationInfo.txt','w')
  simInfoFile.writelines(simInfo)
  simInfoFile.close()


#Function Evaluators
FunctionEval = dgFunctionEvaluator(groups, currentFunction)
BathEval = dgFunctionEvaluator(groups, bathDC.getFunction())
ElevationEval = dgFunctionEvaluator(groups, solution.getFunction())
WindUVEval = dgFunctionEvaluator(groups, windUV)
res = fullMatrixDouble()
#Evaluate Bathymetry Evaluator
outBathEval = open(outputDir+'/Evaluator/evalBathymetry_StdLocations.dat',"w")
outBathGBROOSEval = open(outputDir+'/Evaluator/evalBathymetry_GBROOSLocations.dat',"w")
for point in range(0,len(evalPointsLatLon)):
  #Only write to file if we are on correct cpu:
  BathEval.compute(evalPointsCart[point][0],evalPointsCart[point][1],evalPointsCart[point][2], res)
  outBathEval.write(str(res(0,0))+'\n')
outBathEval.close()
for point in range(0,len(evalPointsGBROOSLatLon)):
  #Only write to file if we are on correct cpu:
  BathEval.compute(evalPointsGBROOSCart[point][0],evalPointsGBROOSCart[point][1],evalPointsGBROOSCart[point][2], res)
  outBathGBROOSEval.write(str(res(0,0))+'\n')
outBathGBROOSEval.close()
# **********************************************
exporterSol.exportIdx(0, t)
# ****************** Iterate  ******************
t_exportStart = Ti + t_exportStart
startcpu = time.clock()

#TEST: OUTPUT VALUE OF WIND (PART 1 OF 2)
if (Msg.GetCommRank() == 0):
  for siteId in range(0,len(evalPointsGBROOSLatLon)):
    outWind = open(outputDir+'/Evaluator/WINDUV_%02d'%(evalPointsGBROOSLatLon[siteId][3])+'.dat',"w")
    outWind.close()

# OPTIONAL 3/3: if starting from an initial solution, adjust start & end steps accordingly
firstStep = 1
lastStep = 1
if (startFromExistingSolution == 1):
  firstStep = stepNbOfInitialSolution + 1
  lastStep = stepNbOfInitialSolution + nbSteps + 1
else:
  firstStep = 1
  lastStep = nbSteps + 1

if (Msg.GetCommRank() == 0):
  print 'Starting simulation from iteration #',firstStep,' ...'
  print ''

for i in range(firstStep, lastStep):
  if(i == lastStep-1):
    dt=Tf-t
  norm = rk.iterate (solution, dt, t)
  t = t +dt
  timeFunction.set(t)
  if ( i%iter == 0 ):
     norm=solution.norm()
     if(Msg.GetCommRank() == 0):
       print''
       print(colored('|ITER| %d of %d',  "red")%(i,lastStep))
       print(colored('------------------------------------------------------', "red"))
       print'|TIME|', printDate(t)
       print'|TIME ELAPSED|', printTime((t - Ti))
       print'| DT |','%.2f'%(dt)
       print'|NORM|','%.4f'%(norm)
       print'|CPUT|',printTime(time.clock() - startcpu)
       print(colored('------------------------------------------------------', "red"))
       print''
  if ( i%export  == 0 and t > t_exportStart ):
    #solution.exportFunctionMsh(currentFunction, outputDir+'/Current/current-%06d'%(i), t, nbExport, '')
    #solution.exportMsh(outputDir+"/Solution/sol-%06d"%(i),t,nbExport)
    exporterSol.exportIdx(i, t)
  if ( i%export_extraTerms == 0 ):
    #exporterDissipationTerm.exportIdx(i, t)
    if(gbrWind!=0):
      windUVDC.interpolate(windUV)
      exporterWindUV.exportIdx(i, t)
      windSourceDC.interpolate(source)
      exporterWindSourceAcc.exportIdx(i, t)
    if(gbrResidual == 2):
      blinkSolDC.interpolate(blink_VelAndElev)
      exporterBlinkSol.exportIdx(i, t)
    if(gbrResidual == 3):
      boundaryGBROOS_MYRSol_DC.interpolate(boundaryUVElev_MYR)
      exporterGBROOS_MYRSol.exportIdx(i, t)
      boundaryGBROOS_CCHSol_DC.interpolate(boundaryUVElev_CCH)
      exporterGBROOS_CCHSol.exportIdx(i, t)
  if ( i%export_fnEval  == 0 ):    
    #Export Function Evaluator results:
    nbExport = nbExport + 1
    #*** Current evaluator: ***
    for point in range(0,len(evalPointsCart)):
      #Only write to file if we are on correct cpu:
      evalPointSP = SPoint3(evalPointsCart[point][0],evalPointsCart[point][1],evalPointsCart[point][2])
      evalPointElem = model.getMeshElementByCoord(evalPointSP)
      whichPartition = evalPointElem.getPartition()
      FunctionEval.compute(evalPointsCart[point][0],evalPointsCart[point][1],evalPointsCart[point][2], res)
      #Project velocities onto lonLat space
      res_uv = fullMatrixDouble(1,2)
      xyVector(res_uv, evalPointsCart[point][0], evalPointsCart[point][1], res)
      if (Msg.GetCommRank() == whichPartition):
        outEval = open(outputDir+'/Evaluator/evalCurrent_%02d'%(point)+'.dat',"a")
        outEval.write(str(t/3600-Ti/3600)+'\t'+str(res.get(0,0))+'\t'+str(res.get(0,1))+'\n')
        outEval.close()
    #*** Elevation evaluator - GBROOS: ***
    for point in range(0,len(evalPointsGBROOSLatLon)):
      #Only write to file if we are on correct cpu:
      evalPointSP = SPoint3(evalPointsGBROOSCart[point][0],evalPointsGBROOSCart[point][1],evalPointsGBROOSCart[point][2])
      evalPointElem = model.getMeshElementByCoord(evalPointSP)
      whichPartition = evalPointElem.getPartition()
      ElevationEval.compute(evalPointsGBROOSCart[point][0],evalPointsGBROOSCart[point][1],evalPointsGBROOSCart[point][2], res)
      res2 = fullMatrixDouble()
      FunctionEval.compute(evalPointsGBROOSCart[point][0],evalPointsGBROOSCart[point][1],evalPointsGBROOSCart[point][2], res2)
      if (Msg.GetCommRank() == whichPartition):
        outElevEval = open(outputDir+'/Evaluator/evalElevGBROOS_%02d'%(evalPointsGBROOSLatLon[point][3])+'.dat',"a")
        outElevEval.write(str(t/3600-Ti/3600)+'\t'+str(res.get(0,0))+'\n')
        outElevEval.close()
        outEval = open(outputDir+'/Evaluator/evalCurrentGBROOS_%02d'%(evalPointsGBROOSLatLon[point][3])+'.dat',"a")
        outEval.write(str(t/3600-Ti/3600)+'\t'+str(res2.get(0,0))+'\t'+str(res2.get(0,1))+'\n')
        outEval.close()
        
      #TEST: OUTPUT VALUE OF WIND (PART 2 OF 2)
      res3 = fullMatrixDouble()
      WindUVEval.compute(evalPointsGBROOSCart[point][0],evalPointsGBROOSCart[point][1],evalPointsGBROOSCart[point][2], res3)
      if (Msg.GetCommRank() == whichPartition):
        outWind = open(outputDir+'/Evaluator/WINDUV_%02d'%(evalPointsGBROOSLatLon[point][3])+'.dat',"a")
        outWind.write(str(t/3600-Ti/3600)+'\t'+str(res3.get(0,0))+'\t'+str(res3.get(0,1))+'\n')
        outWind.close()
        
    #*** Elevation evaluator - Seafarer: ***
    for point in range(0,len(tidalPointsSeafarerLatLon)):
      #Only write to file if we are on correct cpu:
      evalPointSP = SPoint3(tidalPointsSeafarerCart[point][0],tidalPointsSeafarerCart[point][1],tidalPointsSeafarerCart[point][2])
      evalPointElem = model.getMeshElementByCoord(evalPointSP)
      whichPartition = evalPointElem.getPartition()
      ElevationEval.compute(tidalPointsSeafarerCart[point][0],tidalPointsSeafarerCart[point][1],tidalPointsSeafarerCart[point][2], res)
      if (Msg.GetCommRank() == whichPartition):
        outElevEval = open(outputDir+'/Evaluator/evalElevSeafarer_%02d'%(point)+'.dat',"a")
        outElevEval.write(str(t/3600-Ti/3600)+'\t'+str(res.get(0,0))+'\n')
        outElevEval.close()
    #*** Reef Points: current & elevation evaluator: ***
    for point in range(0,len(reefPointsLatLon)):
      #Only write to file if we are on correct cpu:
      evalPointSP = SPoint3(reefPointsCart[point][0],reefPointsCart[point][1],reefPointsCart[point][2])
      evalPointElem = model.getMeshElementByCoord(evalPointSP)
      whichPartition = evalPointElem.getPartition()
      ElevationEval.compute(reefPointsCart[point][0],reefPointsCart[point][1],reefPointsCart[point][2], res)
      if (Msg.GetCommRank() == whichPartition):
        outEval = open(outputDir+'/Evaluator/evalReefPoint_%02d'%(point)+'.dat',"a")
        outEval.write(str(t/3600-Ti/3600)+'\t'+str(res.get(0,0))+'\t'+str(res.get(0,1))+'\t'+str(res.get(0,2))+'\n')
        outEval.close()

endcpu=time.clock()
norm=solution.norm()
Msg.Barrier()
if (Msg.GetCommRank() == 0):
  print ''
  print(colored('********************     End      ********************', "blue"))
  print ''
  print '     - Final Time (Tf):',  printDate(t)
  print '     - Norm of Solution at Tf:',  norm
  print '     - Total CPU Time:', endcpu-startcpu
  print ''
  print(colored('******************************************************', "blue"))
  print ''
Msg.Exit(0)