tutorial_global_oce_in_p/input/gendata.py

0320e25227 Mart*0001 #!/usr/bin/env python
                0002 # -*- coding: iso-8859-15 -*-
                0003 ######################## -*- coding: utf-8 -*-
                0004
                0005 # simple script to generate p-coordinate specific input from standard experiment
                0006 # but unfortunately this script does not reproduce the values in
                0007 # tutorial_global_oce_in_p, yet
                0008
                0009 import numpy as np
                0010 import sys, os
                0011 # requires that the path contains utils/python/MITgcmutils or that the utils
                0012 # are installed via pip or similar:
                0013 import MITgcmutils as mit
                0014
                0015 # some helper routines
                0016 def sq(a):
                0017     a = np.squeeze(a)
                0018     masked_array=np.ma.masked_where(a==0., a)
                0019     return masked_array
                0020
                0021 def readfield(fname,dims,datatype):
                0022     """Call signatures::
                0023
                0024     readfield(filename, dims, numpy.datatype)
                0025
                0026     Read unblocked binary data with dimentions "dims".
                0027     """
                0028
                0029     try:
                0030         fid = open(fname,"rb")
                0031     except:
                0032         sys.exit( fname+": no such file or directory")
                0033
                0034     v   = np.fromfile(fid, datatype)
                0035     fid.close()
                0036
                0037     if sys.byteorder == 'little': v.byteswap(True)
                0038
                0039     if   len(v) == np.prod(dims):     v = v.reshape(dims)
                0040     elif len(v) == np.prod(dims[1:]): v = v.reshape(dims[1:])
                0041     else:
                0042         errstr = (  "dimensions do not match: \n len(data) = " + str(len(v))
                0043                   + ", but prod(dims) = " + str(np.prod(dims)) )
                0044         raise RuntimeError(errstr)
                0045
                0046     return v
                0047
                0048 def writefield(fname,data):
                0049     """Call signatures::
                0050
                0051     writefield(filename, numpy.ndarray)
                0052
                0053     Write unblocked binary data.
                0054     """
                0055
                0056     if True: pass
                0057     else:
                0058         if sys.byteorder == 'little': data.byteswap(True)
                0059
                0060         fid = open(fname,"wb")
                0061         data.tofile(fid)
                0062         fid.close()
                0063
                0064         # switch back to machine format
                0065         if sys.byteorder == 'little': data.byteswap(True)
                0066
                0067 def calc_hydrostatic_pressure(s,t,p0,dz,gravity=9.81,rhoConst=1035.):
                0068     from MITgcmutils import jmd95
                0069     mskz = np.copy(t)
                0070     mskz[mskz!=0]=1.
                0071     dp = np.copy(p0)
                0072     dims=np.asarray(t.shape)
                0073     dims[0]=dims[0]+1
                0074     pf = np.zeros(dims)
                0075     grho = gravity*rhoConst
                0076     rhoInSitu0 = jmd95.dens(s,t,p0/grho)*mskz
                0077     # integration of non-linear hydrostatic equation requires iteration:
                0078     resid = 1
                0079     while resid>1e-15:
                0080         rhoInSitu = jmd95.dens(s,t,p0/grho)*mskz
                0081         # save old pressure
                0082         dp = np.copy(p0)
                0083         # compute new pressure
                0084         pf[0,...] = 0.
                0085         for k in range(nr):
                0086             dpk =  dz[k,...]*gravity*rhoInSitu[k,...]
                0087             p0[k,...]   = (pf[k,...] + 0.5*dpk)*mskz[k,...]
                0088             pf[k+1,...] = (p0[k,...] + 0.5*dpk)*mskz[k,...]
                0089
                0090         # check convergence
                0091         dp = dp-p0
                0092         resid = np.sqrt((dp**2).sum())
                0093         print('hydrostatic pressure: resdiual = %e, '%np.sqrt((dp**2).sum()))
                0094
                0095     print()
                0096     return p0, pf, rhoInSitu
                0097
                0098 gravity = 9.81
                0099 rhoConst= 1035.
                0100 grho=gravity*rhoConst
                0101 nr=15
                0102 ny=40
                0103 nx=90
                0104
                0105 # bathymetry
                0106 prec = 'float64'
                0107 prec = 'float32'
                0108 bdir = '../../tutorial_global_oce_latlon/input'
                0109 bz=np.float64(readfield(os.path.join(bdir,'bathymetry.bin'),[ny,nx],prec))
                0110 # writefield('topog.bin',-(b)*grho)
                0111 b = readfield('topog.bin',[ny,nx],'float64')
                0112 dgp = readfield('deltageopotjmd95.bin',[12,ny,nx],'float64')[0,:,:]
                0113
                0114 # turn fields upside down again for integration in z-levels
                0115 t=readfield('lev_t.bin',[nr,ny,nx],'float64')[::-1,:,:]
                0116 s=readfield('lev_s.bin',[nr,ny,nx],'float64')[::-1,:,:]
                0117 mskz = mit.rdmds('../run/hFacC')[::-1,:,:]
                0118 rac  = mit.rdmds('../run/RAC')
                0119 mskz[mskz!=0] = 1.
                0120 # create geopotential anomaly file:
                0121 delz = np.asarray([50., 70., 100., 140., 190.,
                0122                    240., 290., 340., 390., 440.,
                0123                    490., 540., 590., 640., 690.])
                0124 delp = delz*grho
                0125 # this is from "data", but it is unclear where these numbers come from
                0126 delp = np.asarray([7103300.720021, 6570548.440790, 6041670.010249,
                0127       5516436.666057, 4994602.034410, 4475903.435290,
                0128       3960063.245801, 3446790.312651, 2935781.405664,
                0129       2426722.705046, 1919291.315988, 1413156.804970,
                0130       1008846.750166,  705919.025481,  504089.693499])[::-1]
                0131
                0132 # integrate initial fields vertically minus reference potential
                0133 pf0 = np.hstack([0,np.cumsum(delp)])
                0134 pc0 = 0.5*(pf0[:-1]+pf0[1:])
                0135 pc3d = np.tile(pc0.reshape((nr,1,1)),(1,ny,nx))*mskz
                0136 dp3d = np.tile(delp.reshape((nr,1,1)),(1,ny,nx))
                0137 dz3d = np.tile(delz.reshape((nr,1,1)),(1,ny,nx))
                0138 # first guess of hydrostatic pressure at center-points based on delp
                0139 # pf is the hydrostatic pressure at interfaces (w-points)
                0140 pc = np.copy(pc3d)
                0141 pc,pf,rhoInSitu = calc_hydrostatic_pressure(s,t,pc,dz3d)
                0142
                0143 # the new pressure also implies different delp, here computed as an average
                0144 # over the model domain
                0145 pm = np.zeros((nr,))
                0146 tm = np.zeros((nr,))
                0147 sm = np.zeros((nr,))
                0148 rhom = np.zeros((nr,))
                0149 for k in range(nr):
                0150     racz = rac*mskz[k,:,:]
                0151     pm[k] = (pc[k,:,:]*racz).sum()/racz.sum()
                0152     tm[k] = (t[k,:,:]*racz).sum()/racz.sum()
                0153     sm[k] = (s[k,:,:]*racz).sum()/racz.sum()
                0154     rhom[k] = (rhoInSitu[k,:,:]*racz).sum()/racz.sum()
                0155
                0156 # hydrostatic pressure from averaged temperature and salinity profiles
                0157 pmm = np.copy(pm)
                0158 pmm,pff,rr = calc_hydrostatic_pressure(sm,tm,pmm,delz)
                0159 # this is very similar to diff(pfm), see below
                0160 dp = np.diff(pff)
                0161 print('hydrostatic pressure layer thickness from averaged hydrography:')
                0162 print(' delR = %14f, %14f, %14f,'%(dp[-1],dp[-2],dp[-3]))
                0163 for k in range(nr-4,0,-3):
                0164     print('        %14f, %14f, %14f,'%(dp[k],dp[k-1],dp[k-2]))
                0165
                0166 # averaged pressure at interfaces to compute delP
                0167 pfm = np.zeros((nr+1,))
                0168 for k in range(nr):
                0169     racz = rac*mskz[k,:,:]
                0170     pfm[k] = (pf[k,:,:]*racz).sum()/racz.sum()
                0171
                0172 pfm[nr] = (pf[nr,:,:]*racz).sum()/racz.sum()
                0173 dp = np.diff(pfm)
                0174
                0175 print('hydrostatic pressure layer thickness from averaged pressure:')
                0176 print(' delR = %14f, %14f, %14f,'%(dp[-1],dp[-2],dp[-3]))
                0177 for k in range(nr-4,0,-3):
                0178     print('        %14f, %14f, %14f,'%(dp[k],dp[k-1],dp[k-2]))
                0179
                0180 # now we would like to compute delRc (distance between c-points)
                0181 dp = np.zeros((nr,))
                0182 dp[0]  = pm[0] # assuming zero surface pressure
                0183 dp[1:] = pm[1:]-pm[:-1]
                0184
                0185 print(' delRc = %14f, %14f, %14f,'%(dp[-1],dp[-2],dp[-3]))
                0186 for k in range(nr-4,0,-3):
                0187     print('         %14f, %14f, %14f,'%(dp[k],dp[k-1],dp[k-2]))
                0188
                0189 dp3d = np.tile(dp.reshape((nr,1,1)),(1,ny,nx))
                0190 # this is the correct way of computing the geopotential anomaly
                0191 # (if integr_geoPot = 1)
                0192 geopotanom = -((1./sq(rhoInSitu) - 1/rhoConst)*mskz*dp3d).sum(axis=0)
                0193 # this is equivalent
                0194 geopotanom1= b/rhoConst-(1./sq(rhoInSitu)*mskz*dp3d).sum(axis=0)
                0195 # these are approximation that are not quite accurate
                0196 geopotanom2= ((rhoInSitu - rhoConst)*mskz*dz3d).sum(axis=0)*gravity/rhoConst
                0197 geopotanom3= -((1./sq(rhoInSitu) - 1/rhoConst)*grho*mskz*dz3d).sum(axis=0)
                0198
                0199 # the correct version
                0200 writefield('geopotanom.bin',geopotanom)
                0201
                0202 # plot field
                0203 import matplotlib.pyplot as plt
                0204
                0205 xg = mit.rdmds('../run/XG')
                0206 yg = mit.rdmds('../run/YG')
                0207 f1=plt.figure()
                0208 f1.clf()
                0209 plt.pcolormesh(xg,yg,geopotanom)
                0210 plt.colorbar()
                0211
                0212
                0213 f2=plt.figure()
                0214 f2.clf()
                0215 kindex = -(np.arange(15)+1)
                0216 plt.plot((delp/delz)/gravity,kindex,'x-',label='from delR in "data"')
                0217 plt.plot((np.diff(pff)/delz)/gravity,kindex,':',label='from mean T/S profile')
                0218 plt.plot((np.diff(pfm)/delz)/gravity,kindex,'-.',label='from mean density profile')
                0219 plt.ylabel('k-level')
                0220 plt.title('(delP/delZ)/gravity')
                0221 plt.grid()
                0222 plt.plot([rhoConst,rhoConst],[-1,-15],'--')
                0223 plt.text(rhoConst,-15,'rhoConst=1035.',rotation=270)
                0224 plt.legend()
                0225
                0226 plt.show()