import numpy
from paterson import *
from SetMarineIceSheetBC import *

#Ok, start defining model parameters here

md.timestepping.time_step=0

print "      creating thickness"
h=1000.
md.geometry.thickness=h*numpy.ones((md.mesh.numberofvertices,1))
md.geometry.bed=-1000.*numpy.ones((md.mesh.numberofvertices,1))
md.geometry.surface=md.geometry.bed+md.geometry.thickness;

print "      creating velocities"
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))

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 temperatures"
md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1))

print "      creating flow law parameter"
md.materials.rheology_B=paterson(md.initialization.temperature)
md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))

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

#Deal with boundary conditions:

print "      boundary conditions for diagnostic model"
md=SetMarineIceSheetBC(md,'../Exp/SquareFront.exp')

print "      boundary conditions for thermal model"
md.thermal.spctemperature[:]=md.initialization.temperature
md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1)) 
pos=numpy.nonzero(md.mask.elementongroundedice)[0]
md.basalforcings.geothermalflux[md.mesh.elements[pos,:]-1]=1.*10**-3    #1 mW/m^2
