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