source: issm/trunk/test/Par/ISMIPE.py@ 16560

import numpy
import netCDF4
from SetIceSheetBC import SetIceSheetBC

#Ok, start defining model parameters here

print "      creating thickness"
f = netCDF4.Dataset('../Data/ISMIPE.nc','r')
data = f.variables['data'][:]
f.close()
md.geometry.surface=numpy.zeros((md.mesh.numberofvertices,1))
md.geometry.bed=numpy.zeros((md.mesh.numberofvertices,1))
for i in xrange(0,md.mesh.numberofvertices):
        y=md.mesh.y[i]
        point1=numpy.floor(y/100.)
        point2=numpy.minimum(point1+1,50)
        coeff=(y-(point1-1.)*100.)/100.
        md.geometry.bed[i]=(1.-coeff)*data[point1,1]+coeff*data[point2,1]
        md.geometry.surface[i]=(1.-coeff)*data[point1,2]+coeff*data[point2,2]
md.geometry.thickness=md.geometry.surface-md.geometry.bed
md.geometry.thickness[numpy.nonzero(numpy.logical_not(md.geometry.thickness))]=0.01
md.geometry.bed=md.geometry.surface-md.geometry.thickness

print "      creating drag"
md.friction.coefficient=numpy.zeros((md.mesh.numberofvertices,1))
md.friction.p=numpy.ones((md.mesh.numberofelements,1))
md.friction.q=numpy.ones((md.mesh.numberofelements,1))

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

#Damage
md.damage.D=numpy.zeros((md.mesh.numberofvertices,1))

print "      boundary conditions for stressbalance model:"
#Create node on boundary first (because we can not use mesh)
md=SetIceSheetBC(md)