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))    #q=1.
#Take care of iceshelves: no basal drag
pos=numpy.nonzero(md.mask.elementonfloatingice)[0]
md.friction.coefficient[md.mesh.elements[pos,:]-1]=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.prognostic.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
md.prognostic.spcthickness[pos]=500.
md.prognostic.stabilization=0    #Better result with no artificial diffusivity
md.thermal.stabilization=0
md.timestepping.final_time=500.
md.timestepping.time_step=1
