import numpy
from SetMarineIceSheetBC import *

#Ok, start defining model parameters here

print "      creating thickness"
ymin=numpy.min(md.mesh.y)
ymax=numpy.max(md.mesh.y)
md.geometry.thickness=500.*numpy.ones((md.mesh.numberofvertices,1))
md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness
md.geometry.surface=md.geometry.bed+md.geometry.thickness

print "      creating drag"
md.friction.coefficient=200.*numpy.ones((md.mesh.numberofvertices,1))
md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
md.friction.p=numpy.ones((md.mesh.numberofelements,1))
md.friction.q=numpy.ones((md.mesh.numberofelements,1))

print "      creating initial values"
md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1))
md.initialization.vx=numpy.zeros((md.mesh.numberofvertices,1))
md.initialization.vy=numpy.zeros((md.mesh.numberofvertices,1))
md.initialization.vz=numpy.zeros((md.mesh.numberofvertices,1))
md.initialization.vel=numpy.zeros((md.mesh.numberofvertices,1))
md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1))

print "      creating flow law parameter"
md.materials.rheology_B=1.7687*10**8*numpy.ones((md.mesh.numberofvertices,1))
md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))

print "      creating surface mass balance"
md.surfaceforcings.mass_balance=0.2*numpy.ones((md.mesh.numberofvertices,1))    #0m/a
md.basalforcings.melting_rate=0.*numpy.ones((md.mesh.numberofvertices,1))    #0m/a

print "      boundary conditions"
md=SetMarineIceSheetBC(md,'../Exp/SquareFrontEISMINT.exp')

#Evolution of the ice shelf
pos=numpy.nonzero(md.mesh.y==200000.)    #nodes on the upper boundary condition
md.balancethickness.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
md.balancethickness.spcthickness[pos]=500.
md.masstransport.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
md.masstransport.spcthickness[pos]=500.
md.masstransport.stabilization=0    #Better result with no artificial diffusivity
md.thermal.stabilization=0
md.timestepping.final_time=500.
md.timestepping.time_step=1
