In [ ]:

    

import numpy
import netCDF4
from hashlib import sha1
import matplotlib.pyplot as plt
%matplotlib inline


    

In [ ]:

    

!make convert_Davies_2013/ggge20271-sup-0003-Data_Table1_Eq_lon_lat_Global_HF.nc


    

In [ ]:

    

nc = netCDF4.Dataset('convert_Davies_2013/ggge20271-sup-0003-Data_Table1_Eq_lon_lat_Global_HF.nc','r')


    

In [ ]:

    

# Read dataset cell-center locations and mean heat flow (W/m2)
lat,lon = nc.variables['lat'][:], nc.variables['lon'][:]
mean_HF = nc.variables['mean_HF']
print('Hash of Davies heat flow: ', sha1(mean_HF[:]).hexdigest())


    

In [ ]:

    

plt.pcolormesh(lon,lat,numpy.log10(mean_HF[:])); plt.colorbar(); plt.clim(-2,0);


    

In [ ]:

    

gf = netCDF4.Dataset('ocean_hgrid.nc','r')


    

In [ ]:

    

# Read super grid node locations
xf,yf = gf.variables['x'][:], gf.variables['y'][:]


    

In [ ]:

    

def linterp(lat,lon,Q,yf,xf):
    ni,nj = numpy.shape(lon)[0], numpy.shape(lat)[0]
    xi, yj = (numpy.mod(xf+180.,360))*(ni/360.), (yf+90.)*(nj/180.) # Floating point coords 0<=xi<=ni, 0<=yj<=nj
    wx,wy = numpy.mod(xi+0.5,1.), numpy.mod(yj+0.5,1.) # Weights
    # Calculate indices of model nodes within dataset (takes advantage of uniform resolution of dataset)
    indw, inds = numpy.floor(xi-0.5).astype(int), numpy.floor(yj-0.5).astype(int) # i,j of node to west, south
    inde, indn = numpy.floor(xi+0.5).astype(int), numpy.floor(yj+0.5).astype(int) # i,j of node to east, north
    indw, inde = numpy.mod(indw, ni), numpy.mod(inde, ni) # Periodic in x
    indn = numpy.minimum(indn, nj-1) # Leads to PCM representation for yj>0.5 in northern row
    return (wy*wx*Q[indn,inde] + (1-wy)*(1-wx)*Q[inds,indw]) + (wy*(1-wx)*Q[indn,indw] + (1-wy)*wx*Q[inds,inde])


    

In [ ]:

    

hff = linterp(lat,lon,mean_HF[:],yf,xf) # Heat flux on fine grid nodes
#plt.pcolormesh(xf,yf, numpy.log10(hff)); plt.colorbar();
#hf = hff[1::2,1::2] # Sub-sample
# Integrate with Trapezoidal rule to model grid
hf = 0.5*hff[1::2,:]+0.25*(hff[:-1:2,:]+hff[2::2,:])
hf = 0.5*hf[:,1::2]+0.25*(hf[:,:-1:2]+hf[:,2::2])
print('Hash of re-gridded heat flow: ', sha1(hf).hexdigest())
plt.pcolormesh(xf[::2,::2],yf[::2,::2], numpy.log10(hf)); plt.colorbar(); plt.clim(-2,0);


    

In [ ]:

    

with netCDF4.Dataset('geothermal_davies2013_v1.nc', 'w', format='NETCDF3_64BIT') as newfile:
    newfile.title = 'Geothermal heat flow from Davies, 2013, re-gridded to OM4_025'
    newfile.reference = nc.reference
    newfile.reference_url = nc.reference_url
    newfile.history = numpy.compat.asstr(nc.history[:])+'\nRegridded using Jupyter notebook at https://github.com/NOAA-GFDL/MOM6-examples/blob/dev/master/ice_ocean_SIS2/OM4_025/preprocessing/'
    dj = newfile.createDimension('j', hf.shape[0])
    di = newfile.createDimension('i', hf.shape[1])
    vj = newfile.createVariable('j','f',('j',))
    vj.long_name = 'Grid j-index'
    vi = newfile.createVariable('i','f',('i',))
    vi.long_name = 'Grid i-index'
    vh = newfile.createVariable('geothermal_hf','f',('j','i',))
    vh.long_name = mean_HF.long_name
    vh.standard_name = mean_HF.standard_name
    vh.units = mean_HF.units
    vh.cell_methods = mean_HF.cell_methods
    vh.sha1 = sha1(hf).hexdigest()
    vj[:] = numpy.arange(0.5,hf.shape[0])
    vi[:] = numpy.arange(0.5,hf.shape[1])
    vh[:] = hf[:]
!md5sum geothermal_davies2013_v1.nc


    

In [ ]:

    

jc,jf=44,1007
plt.figure(figsize=(12,6))
plt.plot(lon,( (mean_HF[jc,:]+mean_HF[jc+1,:])/2 ) );
plt.plot(numpy.mod( xf[jf,:]+180 ,360)-180,hff[jf,:],'.')
plt.plot(numpy.mod( xf[jf,1::2]+180 ,360)-180,( hf[int(jf/2),:] ) ,'.');