[14134] | 1 | import numpy
|
---|
| 2 | import netCDF4
|
---|
[16560] | 3 | from SetIceSheetBC import SetIceSheetBC
|
---|
[14134] | 4 |
|
---|
| 5 | #Ok, start defining model parameters here
|
---|
| 6 |
|
---|
| 7 | print " creating thickness"
|
---|
| 8 | f = netCDF4.Dataset('../Data/ISMIPE.nc','r')
|
---|
| 9 | data = f.variables['data'][:]
|
---|
| 10 | f.close()
|
---|
| 11 | md.geometry.surface=numpy.zeros((md.mesh.numberofvertices,1))
|
---|
[17806] | 12 | md.geometry.base=numpy.zeros((md.mesh.numberofvertices,1))
|
---|
[14134] | 13 | for i in xrange(0,md.mesh.numberofvertices):
|
---|
| 14 | y=md.mesh.y[i]
|
---|
| 15 | point1=numpy.floor(y/100.)
|
---|
| 16 | point2=numpy.minimum(point1+1,50)
|
---|
| 17 | coeff=(y-(point1-1.)*100.)/100.
|
---|
[17806] | 18 | md.geometry.base[i]=(1.-coeff)*data[point1,1]+coeff*data[point2,1]
|
---|
[14134] | 19 | md.geometry.surface[i]=(1.-coeff)*data[point1,2]+coeff*data[point2,2]
|
---|
[17806] | 20 | md.geometry.thickness=md.geometry.surface-md.geometry.base
|
---|
[14134] | 21 | md.geometry.thickness[numpy.nonzero(numpy.logical_not(md.geometry.thickness))]=0.01
|
---|
[17806] | 22 | md.geometry.base=md.geometry.surface-md.geometry.thickness
|
---|
[14134] | 23 |
|
---|
| 24 | print " creating drag"
|
---|
| 25 | md.friction.coefficient=numpy.zeros((md.mesh.numberofvertices,1))
|
---|
| 26 | md.friction.p=numpy.ones((md.mesh.numberofelements,1))
|
---|
| 27 | md.friction.q=numpy.ones((md.mesh.numberofelements,1))
|
---|
| 28 |
|
---|
| 29 | print " creating flow law parameter"
|
---|
| 30 | md.materials.rheology_B=6.8067*10**7*numpy.ones((md.mesh.numberofvertices,1))
|
---|
| 31 | md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
|
---|
| 32 |
|
---|
[16137] | 33 | print " boundary conditions for stressbalance model:"
|
---|
[14134] | 34 | #Create node on boundary first (because we can not use mesh)
|
---|
| 35 | md=SetIceSheetBC(md)
|
---|