isomip/input.obcs/gendata.py

abfe198bce Mart*0001 # simple python script to generate additional input data
                0002 import numpy as np
                0003 import sys
                0004
                0005 kwr, kprt =1, 0
                0006 nx, ny, nr = 50, 100, 30
                0007
                0008 def writefield(fname,data):
                0009     """Call signatures::
                0010
                0011     writefield(filename, numpy.ndarray)
f18a893d42 Mart*0012
abfe198bce Mart*0013     Write unblocked binary data.
                0014     """
                0015
                0016     if sys.byteorder == 'little': data.byteswap(True)
                0017
                0018     fid = open(fname,"wb")
                0019     data.tofile(fid)
                0020     fid.close()
                0021
                0022     # go back to original byte ordering
                0023     if sys.byteorder == 'little': data.byteswap(True)
                0024
                0025 def readfield(fname,dims,datatype):
                0026     """Call signatures::
                0027
                0028     readfield(filename, dims, numpy.datatype)
f18a893d42 Mart*0029
abfe198bce Mart*0030     Read unblocked binary data with dimentions "dims".
                0031     """
                0032     fid = open(fname,"rb")
                0033     v   = np.fromfile(fid, datatype)
                0034     fid.close()
                0035
                0036     if sys.byteorder == 'little': v.byteswap(True)
                0037
                0038     if   len(v) == np.prod(dims):     v = v.reshape(dims)
                0039     elif len(v) == np.prod(dims[1:]): v = v.reshape(dims[1:])
                0040     else:
f18a893d42 Mart*0041         errstr = (  "dimensions do not match: \n len(data) = " + str(len(v))
abfe198bce Mart*0042                   + ", but prod(dims) = " + str(np.prod(dims)) )
                0043         raise RuntimeError(errstr)
                0044
                0045     return v
                0046
                0047 y0 = 50
                0048 # use default topography and modify it
                0049 fld = readfield("../input/icetopo.exp1",[ny,nx],'float64')
                0050 topo = np.copy(fld)
                0051 topo[y0:,:] = 0.
                0052 writefield("icetopo.obcs",topo)
                0053
                0054 fld = readfield("../input/iceShelf_Mass.bin",[ny,nx],'float64')
                0055 # new ice mass
                0056 iceMass = np.copy(fld)
                0057 iceMass[y0:,:] = 0.
                0058 writefield("iceShelf_Mass.obcs",iceMass)
                0059 # ice mass tendency
                0060 iceMassTend = np.zeros(iceMass.shape)
                0061 dt = 1800.
                0062 iceMassTend[:40,:]=iceMass[:40,:]*1.e-5/dt
                0063 writefield("iceShelf_MassTend.obcs",iceMassTend)
                0064
                0065 # for averaging
                0066 msk = np.ones(iceMass.shape)
                0067 msk[y0:,:] = 2.
                0068 writefield("iceMask.obcs",msk)
f18a893d42 Mart*0069
                0070 # salinity profile for stability
                0071 dz = 30.
                0072 N2 = 1.e-5
                0073 sBeta = 7e-4
                0074 gravity=9.81
                0075 dSdz = N2/gravity/sBeta * np.ones((30,))
                0076 S0 = 34.4
                0077 sRef0 = np.cumsum(dSdz*dz)
                0078 sRef = np.maximum(sRef0-sRef0[9],0.) + S0
                0079 #sRef = np.maximum(sRef0-sRef0[0],0.) + S0