interp

interp/get_alltypes_cdo_griddes.sh

@@ -0,0 +1,65 @@
+#!/bin/bash
+set -u
+
+# This script generates CDO grid file descriptions.
+# These CDO description files are used for horizontal interpolation from/to a given grid.
+# The script assumes that the horizontal dimensions are called (y,x) but feel free to edit.
+
+# Basically, from the netCDF gridfile, this scripts creates a dummy netCDF that CDO recognizes, to get the CDO grid description file. We'll do that:
+# - longitude variable with netCDF attributes standard_name="longitude", units="degrees East"
+# - latitude variable with netCDF attributes standard_name="latitude", units="degrees North"
+# - one dummy variable with a netCDF attribute coordinates='<lat> <lon>'
+
+# Then cdo griddes dummy.nc > my_grid.grid will work.
+
+# Input file with lat, lon and depth
+input_file="/storepelican/pelletie/nemo_bdy/grids/coord_paraso.nc"
+output_folder="/storepelican/pelletie/nemo_bdy/grids"
+
+module purge
+module load 2018b CDO NCO
+
+pid="$$"
+
+types=('t' 'u' 'v')
+
+mkdir -p ${output_folder}
+cd ${output_folder}
+
+for type in "${types[@]}"; do
+
+    lat_var="gphi${type}"
+    lon_var="glam${type}"
+
+    output_griddes="eORCA025-SO_grid_${type^}.grid"
+    
+    # Just take lat/lon from file
+    ncks -C -v ${lat_var},${lon_var} ${input_file} -O tmp_pid${pid}.nc
+    # Edit lat attributes
+    ncatted -a standard_name,${lat_var},o,c,"latitude" -a units,${lat_var},o,c,"degrees North" tmp_pid${pid}.nc
+    # Edit lon attributes
+    ncatted -a standard_name,${lon_var},o,c,"longitude" -a units,${lon_var},o,c,"degrees East" tmp_pid${pid}.nc
+    # Create some dummy variable. Here, dimension names are assumed to be y/x. Change if needed.
+    ncap2 -O -s "dummy[\$y,\$x]=1.f;" tmp_pid${pid}.nc tmp2_pid${pid}.nc
+    # Edit dummy variable attributes to say its coordinates are you lat/lons.
+    ncatted -a coordinates,dummy,o,c,"${lat_var} ${lon_var}" tmp2_pid${pid}.nc
+
+    # Get CDO griddes
+    cdo griddes tmp2_pid${pid}.nc > ${output_griddes}
+
+done
+
+rm -f ./*pid${pid}*
+
+
+
+# Expert level: if you want to perform conservative interpolation, you'll need grid corners as well.
+# For ORCA grids (the NEMO grids), this can be done using F nodes (grid corners), check `get_orca_corners.py`
+# Then, you need:
+# 1) One additional dimension <corn> with size 4 (the name does not matter) and two additional fields:
+#   1a) <lon_corn>(y,x,<corn>): longitudes of grid corners
+#   1b)  <lat_corn>(y,x,<corn>): latitudes of grid corners
+# 2) Add an attributes "bounds" to the lat/lon fields in the dummy file, i.e.:
+#     ncatted -a standard_name,${lat_var},o,c,"latitude" -a units,${lat_var},o,c,"degrees North" -a bounds,${lat_var},o,c,<lat_corn> tmp_pid${pid}.nc
+#    ncatted -a standard_name,${lon_var},o,c,"longitude" -a units,${lon_var},o,c,"degrees East" -a bounds,${lon_var},o,c,<lon_corn> tmp_pid${pid}.nc
+# Then, create the dummy var the same way as above. The CDO gridfile description should then be much larger (containing corners). And can be used for conservative interpolation.

interp/get_bisicles_corners.py

@@ -0,0 +1,122 @@
+import sys
+import numpy as np
+from netCDF4 import Dataset
+import os
+import time
+from pyproj import CRS, Transformer
+
+# C. Pelletier july 2021.
+
+# Example of using PYPROJ to manage projections.
+# I highly recommend using this tool, as you're then relying on something documented from projecting (i.e. from lat/lon to stereographic)
+# This is pretty useful for ice sheet / ocean coupling and pre/post processing. since most ice sheet models are defined on polar stereographic grids.
+
+# This script just gets the grid corners of the BISICLE grid.
+# But it's an illustration of how to use PYPROJ.
+
+# <https://pyproj4.github.io/pyproj/stable/>
+# <https://proj.org/>
+
+grid_bisicles = "/afast/pelletie/totten/BISICLES-smaller_grid.nc"
+output_file = "/afast/pelletie/totten/BISICLES-smaller_grid_corners.nc"
+
+# description of the projection.
+proj_polster = "+proj=stere +lat_0=-90.0 +lat_ts=-71. +lon_0=-90."
+# The projection python object
+crs_polster = CRS.from_proj4(proj_polster)
+
+tmpdat = Dataset(grid_bisicles, mode='r')
+
+# Center of cells in stereographic projection. No idea why X and Y are interverted.
+x_ctn = tmpdat.variables['Y'][:]
+y_ctn = tmpdat.variables['X'][:]
+
+mask = tmpdat.variables['MASK'][:].astype(int)
+
+tmpdat.close()
+
+nx = np.size(x_ctn)
+ny = np.size(y_ctn)
+
+
+# Edges of cell in stereographics projection. This is just the middle of two consecutive cells. And a bit of symmetry for the edges.
+x_stg = np.ndarray(shape=[nx+1], dtype=float)
+y_stg = np.ndarray(shape=[ny+1], dtype=float)
+
+x_stg[1:-1] = .5 * ( x_ctn[:-1] + x_ctn[1:] )
+x_stg[0] = 2 * x_ctn[0] - x_stg[1]
+x_stg[-1] = 2 * x_ctn[-1] - x_stg[-2]
+
+y_stg[1:-1] = .5 * ( y_ctn[:-1] + y_ctn[1:] )
+y_stg[0] = 2 * y_ctn[0] - y_stg[1]
+y_stg[-1] = 2 * y_ctn[-1] - y_stg[-2]
+
+
+# PYPROJ transformation object, polster -> latlon
+# <crs_something>.geodetic_crs is always the latlon projection, regardless of <crs_something>
+polster_to_latlon = Transformer.from_crs(crs_polster, crs_polster.geodetic_crs)
+
+# Get lon, lat of centers from x,y (which are broadcasted with [None,:] etc.)
+
+lon_ctn, lat_ctn = polster_to_latlon.transform(x_ctn[None,:] * (np.ones(shape=[ny], dtype=float))[:,None], \
+                                               y_ctn[:,None] * (np.ones(shape=[nx], dtype=float))[None,:])
+
+# No idea why I'm doing this. I guess everything's corked and dirty.
+lon_ctn = - lon_ctn
+
+# Reformat corners in (y,x,4) instead of (y+1,x+1).
+x_crn = np.ndarray(shape=[ny, nx, 4], dtype=float)
+y_crn = np.ndarray(shape=[ny, nx, 4], dtype=float)
+
+x_crn[:,:,0] = (x_stg[1:])[None,:] * (np.ones(shape=[ny], dtype=float))[:,None]
+y_crn[:,:,0] = (y_stg[1:])[:,None] * (np.ones(shape=[nx], dtype=float))[None,:]
+
+x_crn[:,:,1] = (x_stg[1:])[None,:] * (np.ones(shape=[ny], dtype=float))[:,None]
+y_crn[:,:,1] = (y_stg[:-1])[:,None] * (np.ones(shape=[nx], dtype=float))[None,:]
+
+x_crn[:,:,2] = (x_stg[:-1])[None,:] * (np.ones(shape=[ny], dtype=float))[:,None]
+y_crn[:,:,2] = (y_stg[:-1])[:,None] * (np.ones(shape=[nx], dtype=float))[None,:]
+
+x_crn[:,:,3] = (x_stg[:-1])[None,:] * (np.ones(shape=[ny], dtype=float))[:,None]
+y_crn[:,:,3] = (y_stg[1:])[:,None] * (np.ones(shape=[nx], dtype=float))[None,:]
+
+
+
+lon_crn = np.ndarray(shape=[ny, nx, 4], dtype=float)
+lat_crn = np.ndarray(shape=[ny, nx, 4], dtype=float)
+
+for m in range(0,4):
+    lon_crn[:,:,m], lat_crn[:,:,m] = polster_to_latlon.transform(x_crn[:,:,m], y_crn[:,:,m])
+
+
+lon_crn = - lon_crn
+
+out_nc = Dataset(output_file, mode='w')
+
+out_nc.createDimension('y', ny)
+out_nc.createDimension('x', nx)
+out_nc.createDimension('crn', 4)
+
+tmp = out_nc.createVariable('lat', datatype='f8', dimensions=['y', 'x'])
+tmp[:] = lat_ctn
+tmp.standard_name = 'latitude'
+tmp.units = 'degrees North'
+tmp.bounds = 'lat_corn'
+
+tmp = out_nc.createVariable('lat_corn', datatype='f8', dimensions=['y', 'x','crn'])
+tmp[:] = lat_crn
+
+tmp = out_nc.createVariable('lon', datatype='f8', dimensions=['y', 'x'])
+tmp[:] = lon_ctn
+tmp.standard_name = 'longitude'
+tmp.units = 'degrees East'
+tmp.bounds = 'lon_corn'
+
+tmp = out_nc.createVariable('lon_corn', datatype='f8', dimensions=['y', 'x','crn'])
+tmp[:] = lon_crn
+
+tmp = out_nc.createVariable('mask', datatype='i2', dimensions=['y', 'x'])
+tmp[:] = mask
+tmp.coordinates = 'lat lon'
+
+out_nc.close()

interp/get_cdo_griddes.sh

@@ -0,0 +1,55 @@
+#!/bin/bash
+set -u
+
+# This script generates CDO grid file descriptions.
+# These CDO description files are used for horizontal interpolation from/to a given grid.
+# The script assumes that the horizontal dimensions are called (y,x) but feel free to edit.
+
+# Basically, from the netCDF gridfile, this scripts creates a dummy netCDF that CDO recognizes, to get the CDO grid description file. We'll do that:
+# - longitude variable with netCDF attributes standard_name="longitude", units="degrees East"
+# - latitude variable with netCDF attributes standard_name="latitude", units="degrees North"
+# - one dummy variable with a netCDF attribute coordinates='<lat> <lon>'
+
+# Then cdo griddes dummy.nc > my_grid.grid will work.
+
+# Input file with lat, lon and depth
+input_file="/mfast/pelletie/oras5/grids/mesh_mask.nc"
+lat_var="gphit"
+lon_var="glamt"
+depth_var="gdept_1d"
+
+output_griddes="/mfast/pelletie/oras5/grids/example_mine.grid"
+
+module purge
+module load 2018b CDO NCO
+
+pid="$$"
+
+# Just take lat/lon from file
+ncks -C -v ${lat_var},${lon_var} ${input_file} -O tmp_pid${pid}.nc
+# Edit lat attributes
+ncatted -a standard_name,${lat_var},o,c,"latitude" -a units,${lat_var},o,c,"degrees North" tmp_pid${pid}.nc
+# Edit lon attributes
+ncatted -a standard_name,${lon_var},o,c,"longitude" -a units,${lon_var},o,c,"degrees East" tmp_pid${pid}.nc
+# Create some dummy variable. Here, dimension names are assumed to be y/x. Change if needed.
+ncap2 -O -s "dummy[\$y,\$x]=1.f;" tmp_pid${pid}.nc tmp2_pid${pid}.nc
+# Edit dummy variable attributes to say its coordinates are you lat/lons.
+ncatted -a coordinates,dummy,o,c,"${lat_var} ${lon_var}" tmp2_pid${pid}.nc
+
+# Get CDO griddes
+cdo griddes tmp2_pid${pid}.nc > ${output_griddes}
+
+rm -f ./*pid${pid}*
+
+
+
+# Expert level: if you want to perform conservative interpolation, you'll need grid corners as well.
+# For ORCA grids (the NEMO grids), this can be done using F nodes (grid corners), check `get_orca_corners.py`
+# Then, you need:
+# 1) One additional dimension <corn> with size 4 (the name does not matter) and two additional fields:
+#   1a) <lon_corn>(y,x,<corn>): longitudes of grid corners
+#   1b)  <lat_corn>(y,x,<corn>): latitudes of grid corners
+# 2) Add an attributes "bounds" to the lat/lon fields in the dummy file, i.e.:
+#     ncatted -a standard_name,${lat_var},o,c,"latitude" -a units,${lat_var},o,c,"degrees North" -a bounds,${lat_var},o,c,<lat_corn> tmp_pid${pid}.nc
+#    ncatted -a standard_name,${lon_var},o,c,"longitude" -a units,${lon_var},o,c,"degrees East" -a bounds,${lon_var},o,c,<lon_corn> tmp_pid${pid}.nc
+# Then, create the dummy var the same way as above. The CDO gridfile description should then be much larger (containing corners). And can be used for conservative interpolation.

interp/get_orca_corners.py

@@ -0,0 +1,116 @@
+import sys
+
+import numpy as np
+from netCDF4 import Dataset
+
+import os
+import time
+
+
+org_grid = '/afast/pelletie/orca025/eORCA025-SO_BMbathy_surf.nc'
+dest_grid = '/afast/pelletie/orca025/eORCA025-SO_corners_new.nc'
+
+data_org = Dataset(org_grid, mode='r')
+
+glamt = data_org.variables['glamt'][:]
+gphit = data_org.variables['gphit'][:]
+
+glamu = data_org.variables['glamu'][:]
+gphiu = data_org.variables['gphiu'][:]
+
+glamv = data_org.variables['glamv'][:]
+gphiv = data_org.variables['gphiv'][:]
+
+glamf = data_org.variables['glamf'][:]
+gphif = data_org.variables['gphif'][:]
+
+mask = data_org.variables['tmask'][:]
+
+data_org.close()
+
+ny, nx = np.shape(glamt)
+
+out_nc = Dataset(dest_grid, mode='w')
+
+out_nc.createDimension('y', ny)
+out_nc.createDimension('x', nx)
+
+tmp = out_nc.createVariable(varname='lon', datatype='f8', dimensions=['y', 'x'])
+tmp[:] = glamt[:]
+tmp.standard_name = "longitude"
+tmp.units = "degrees East"
+tmp.bounds = "lon_corn"
+
+tmp = out_nc.createVariable(varname='lat', datatype='f8', dimensions=['y', 'x'])
+tmp[:] = gphit[:]
+tmp.standard_name = "latitude"
+tmp.units = "degrees North"
+tmp.bounds = "lat_corn"
+
+tmp = out_nc.createVariable(varname='mask', datatype='b', dimensions=['y', 'x'])
+tmp[:] = mask[:]
+tmp.coordinates = "lat lon"
+
+
+out_nc.createDimension('corners', 4)
+
+corn_tmp = np.ndarray(shape=[ny+1,nx+1], dtype=float)
+
+corn_tmp[1:, 1:] = glamf
+corn_tmp[0,1:] = 2*glamu[0,:] - glamf[0,:]
+corn_tmp[:,0] = corn_tmp[:,-3] #- 360.
+
+lon_c = out_nc.createVariable(varname='lon_corn', datatype='f8', dimensions=['y', 'x', 'corners'])
+
+lon_c[:,:,0] = corn_tmp[:-1,:-1]
+lon_c[:,:,1] = corn_tmp[1:,:-1]
+lon_c[:,:,2] = corn_tmp[:-1,1:]
+lon_c[:,:,3] = corn_tmp[1:,1:]
+
+res = .25
+
+corr_lonc = np.logical_or(np.abs(lon_c[:,:,0] - lon_c[:,:,2]) < res/100.,
+                          np.abs(lon_c[:,:,1] - lon_c[:,:,3]) < res/100.)
+
+lon_c[:,:,0] = np.where(corr_lonc,\
+                        glamt[:] - .5 * res, lon_c[:,:,0])
+lon_c[:,:,1] = np.where(corr_lonc,\
+                        glamt[:] - .5 * res, lon_c[:,:,1])
+lon_c[:,:,2] = np.where(corr_lonc,\
+                        glamt[:] + .5 * res, lon_c[:,:,2])
+lon_c[:,:,3] = np.where(corr_lonc,\
+                        glamt[:] + .5 * res, lon_c[:,:,3])
+
+
+del lon_c
+
+corn_tmp[1:, 1:] = gphif
+corn_tmp[0,1:] = 2*gphiu[0,:] - gphif[0,:]
+corn_tmp[:,0] = corn_tmp[:,-1]
+
+lat_c = out_nc.createVariable(varname='lat_corn', datatype='f8', dimensions=['y', 'x', 'corners'])
+
+lat_c[:,:,0] = corn_tmp[:-1,:-1]
+lat_c[:,:,1] = corn_tmp[1:,:-1]
+lat_c[:,:,2] = corn_tmp[:-1,1:]
+lat_c[:,:,3] = corn_tmp[1:,1:]
+
+corr_latc = np.logical_or(np.abs(lat_c[:,:,0] - lat_c[:,:,1]) < res/100.,
+                          np.abs(lat_c[:,:,2] - lat_c[:,:,3]) < res/100.)
+
+lat_c[:,:,0] = np.where(corr_latc,\
+                        gphit[:] - .5 * res, lat_c[:,:,0])
+lat_c[:,:,1] = np.where(corr_latc,\
+                        gphit[:] + .5 * res, lat_c[:,:,1])
+lat_c[:,:,2] = np.where(corr_latc,\
+                        gphit[:] - .5 * res, lat_c[:,:,2])
+lat_c[:,:,3] = np.where(corr_latc,\
+                        gphit[:] + .5 * res, lat_c[:,:,3])
+
+
+
+out_nc.close()
+
+
+
+

interp/get_temp_sublayers_lucy.py

@@ -0,0 +1,160 @@
+import numpy as np
+from netCDF4 import Dataset
+import os
+import sys
+
+# path to input 4D file    
+input_file = "/afast/pelletie/nemo_outputs/so025_noisf/pped/so025_noisf_1m_tempsal_1999.nc"
+# path to output file
+out_file = "/afast/pelletie/nemo_outputs/so025_noisf/pped/so025_noisf_1m_tempsal_layers_1999.nc"
+
+# Upper and lower bounds of average ocean temperature layers.
+upper_bounds = [0.,   50., 150., 300.]
+lower_bounds = [50., 150., 300., 500.]
+
+# path to 'mesh mask' file
+mesh_file = "/afast/pelletie/nemo_outputs/so025_noisf/mesh/mesh_mask.nc"
+
+
+# end input
+
+tmpdat = Dataset(mesh_file, mode='r')
+e3t = tmpdat.variables['e3t_0'][:]
+gdepth_3d = tmpdat.variables['gdept_0'][:]
+tmask = (tmpdat.variables['tmask'][:]).astype(bool)
+tmpdat.close()
+
+
+nz, ny, nx = np.shape(tmask)
+
+tmask = np.where((np.logical_not(tmask[0,:,:]))[None,:,:] * (np.ones(shape=[nz], dtype=bool))[:,None,None], \
+                 False, tmask)
+
+
+n_layers = np.size(upper_bounds)
+thick_intg = np.nan * np.ndarray(shape=[n_layers,nz, ny, nx], dtype=float)
+# thick_intg = np.zeros(shape=[n_layers,nz, ny, nx], dtype=float)
+
+
+upper_mesh = gdepth_3d - .5 * e3t
+lower_mesh = gdepth_3d + .5 * e3t
+
+for m in range(0,n_layers):
+
+    fillidx = np.nonzero(np.logical_and( upper_bounds[m] <= upper_mesh, lower_bounds[m] >= lower_mesh))
+    (thick_intg[m,:,:,:])[fillidx] = e3t[fillidx]
+    
+    fillidx = np.nonzero(np.logical_or( upper_bounds[m] > lower_mesh, lower_bounds[m] < upper_mesh))
+    (thick_intg[m,:,:,:])[fillidx] = 0.
+
+    fillidx = np.nonzero(np.logical_and( upper_bounds[m] > upper_mesh, upper_bounds[m] <= lower_mesh))
+    (thick_intg[m,:,:,:])[fillidx] = (lower_mesh - upper_bounds[m])[fillidx]
+
+    fillidx = np.nonzero(np.logical_and( lower_bounds[m] < lower_mesh, lower_bounds[m] >= upper_mesh))
+    (thick_intg[m,:,:,:])[fillidx] = (lower_bounds[m] - upper_mesh)[fillidx]
+
+    fillidx = np.nonzero(np.logical_not(tmask))
+    (thick_intg[m,:,:,:])[fillidx] = 0.
+
+    # print('nan count', np.count_nonzero(np.isnan(thick_intg[m,:,:,:])), nz*ny*nx)
+
+
+#     print(m,thick_intg[m,:,300,500], upper_mesh[:,300,500], lower_mesh[:,300,500], \
+#           tmask[:,300,500])
+
+# sys.exit()
+
+    
+
+    
+del upper_mesh,lower_mesh,e3t, gdepth_3d
+
+
+fillval = -1.e20
+
+
+
+tmpdat = Dataset(input_file, mode='r')
+
+in_temp = tmpdat.variables['thetao'][:]
+in_sal  = tmpdat.variables['so'][:]
+
+lat = tmpdat.variables['nav_lat'][:]
+lon = tmpdat.variables['nav_lon'][:]
+depth_bounds = tmpdat.variables['deptht_bounds'][:]
+depth = tmpdat.variables['deptht'][:]
+
+out_nc = Dataset(out_file, mode='w')
+
+out_nc.createDimension('y', ny)
+out_nc.createDimension('x', nx)
+out_nc.createDimension('depth', n_layers)
+out_nc.createDimension('time_counter', 0)
+out_nc.createDimension('axis_nbounds', 2)
+
+tmp = out_nc.createVariable(varname='time_counter', datatype='f8', dimensions=['time_counter'])
+tmp[:] = tmpdat.variables['time_counter'][:]
+nt = np.size(tmp[:])
+tmp.setncatts( tmpdat.variables['time_counter'].__dict__)
+
+tmp = out_nc.createVariable(varname='time_counter_bounds', datatype='f8', dimensions=['time_counter', 'axis_nbounds'])
+tmp[:] = tmpdat.variables['time_counter_bounds'][:]
+tmp.setncatts( tmpdat.variables['time_counter_bounds'].__dict__)
+
+tmp = out_nc.createVariable(varname='nav_lat', datatype='f8', dimensions=['y', 'x'])
+tmp[:] = tmpdat.variables['nav_lat'][:]
+tmp.setncatts( tmpdat.variables['nav_lat'].__dict__)
+
+tmp = out_nc.createVariable(varname='nav_lon', datatype='f8', dimensions=['y', 'x'])
+tmp[:] = tmpdat.variables['nav_lon'][:]
+tmp.setncatts( tmpdat.variables['nav_lon'].__dict__)
+
+tmp = out_nc.createVariable(varname='depth', datatype='f8', dimensions=['depth'])
+tmp[:] = np.mean(np.vstack((upper_bounds, lower_bounds)), axis=0)
+tmp.standard_name = "depth"
+tmp.units="m"
+tmp.axis="Z"
+tmp.positive="down"
+tmp.bounds="depth_bounds"
+
+tmp = out_nc.createVariable(varname='depth_bounds', datatype='f8', dimensions=['depth', 'axis_nbounds'])
+tmp[:,0] = upper_bounds
+tmp[:,1] = lower_bounds
+
+out_temp = out_nc.createVariable(varname='thetao', datatype='f4', dimensions=['time_counter', 'depth', 'y', 'x'], fill_value=fillval)
+out_temp.standard_name = "conservative_temperature_mean_on_sublayers"
+out_temp.units = "degC"
+out_temp.coordinates = "time_counter depth nav_lat nav_lon"
+
+out_sal = out_nc.createVariable(varname='so', datatype='f4', dimensions=['time_counter', 'depth', 'y', 'x'], fill_value=fillval)
+out_sal.standard_name = "absolute_salinity_mean_on_sublayers"
+out_sal.units = "g/kg"
+out_sal.coordinates = "time_counter depth nav_lat nav_lon"
+
+
+
+for m in range(0,n_layers):
+
+    out_temp[:,m,:,:] = np.where( (np.all(thick_intg[m,:,:,:] == 0., axis=0))[None,:,:] * (np.ones(shape=[nt], dtype=bool))[:,None,None], \
+                                 fillval, \
+                                 np.sum(in_temp *   (thick_intg[m,:,:,:])[None,:,:,:] * (np.ones(shape=[nt], dtype=float))[:,None,None,None], axis=1) / np.sum((thick_intg[m,:,:,:])[None,:,:,:] * (np.ones(shape=[nt], dtype=float))[:,None,None,None], axis=1))
+
+
+    out_sal[:,m,:,:] = np.where( (np.all(thick_intg[m,:,:,:] == 0., axis=0))[None,:,:] * (np.ones(shape=[nt], dtype=bool))[:,None,None], \
+                                 fillval, \
+                                 np.sum(in_sal *   (thick_intg[m,:,:,:])[None,:,:,:] * (np.ones(shape=[nt], dtype=float))[:,None,None,None], axis=1) / np.sum((thick_intg[m,:,:,:])[None,:,:,:] * (np.ones(shape=[nt], dtype=float))[:,None,None,None], axis=1))
+
+    
+    # print('dbg thick_intg', m, thick_intg[m,:,301,490])
+    # print('dbg out temp', m, out_temp[m,:,301,490])
+
+    # out_temp[:,m,:,:] = np.sum( (thick_intg[m,:,:,:])[None,:,:,:] * in_temp, axis=1)
+    # print('isnan out', np.count_nonzero(np.isnan(out_temp[0,m,:,:])))
+
+
+# out_nc.createDimension('z', nz)
+# tmp = out_nc.createVariable(varname='thick_intg', datatype='f8', dimensions=['layer', 'z', 'y', 'x'])
+# tmp[:] = thick_intg[:]
+
+tmpdat.close()
+out_nc.close()

interp/nn_interp.py

@@ -0,0 +1,28 @@
+import numpy as np
+
+def nn_interp(input_arr):
+
+    shape = np.shape(input_arr)
+    n_shape = np.size(shape)
+    
+    if(n_shape < 2):
+        print('error. cant do nn_interp on arrays with rank < 2')
+        sys.exit()
+    
+    else:
+
+        iter_arr = np.ma.copy(input_arr)
+        
+        while( np.count_nonzero(iter_arr.mask) > 0 ):
+            for shift in (-1,1):
+                for axis in (-1,-2):        
+                    arr_shifted=np.roll(iter_arr,shift=shift,axis=axis)
+                    idx=~arr_shifted.mask * iter_arr.mask
+                    iter_arr[idx]=arr_shifted[idx]
+
+    return iter_arr
+
+        
+
+
+        

interp/readme.md

@@ -0,0 +1,59 @@
+# CDO horizontal interpolation.
+
+1. Read the [CDO documentation](https://code.mpimet.mpg.de/projects/cdo/).
+
+2. Prepare CDO grid description files (see `get_cdo_griddes.sh`), which are kind of "blueprints" for CDO - it describes your grid.
+
+3. Your input data should be _properly masked_. This means that the masked cells should be `_FillValue`, and appear in white when you look at them with `ncview`. This way, CDO will know that the fields are masked there, and that the values should not be used. If there's one "dummy" value that is _not_ the specified `_FillValue`, then CDO will use these values as if they were good. Which you don't want to do. To properly mask a field, try to do something like that:
+ ```
+ ncatted -a _FillValue,<some_var>,o,f,1.e20 file.nc # in <some_var> is a SIMPLE precision FLOAT
+ ncatted -a _FillValue,<some_var>,o,d,1.e20 file.nc # in <some_var> is a DOUBLE precision float
+ ```
+ And then:
+ ```
+ ncap2 -O -s 'where(mask==0) <some_var>=<some_var>@_FillValue;' file.nc output.nc
+ ```
+ This assumes that `file.nc` has a field called `mask`, and that cells where `mask` is zero ought to be masked out.  `<some_var>` can have more dimensions than `mask`: for example, if `mask(y,x)` and `<some_var>(t,y,x)`, things will be broadcasted along the `t` dimension (all `t` indices will be masked). Think about what you do: if you're masking a 3D `(z,y,x)` ocean field, use a 3D mask; the mask is not the same at the surface and at 5000m depth.
+ Other things can be used in the `where` statement, e.g. `where(abs(<some_var>) > 1.e10)` will mask anywhere where `<some_var>` exceeds `1.e10`.
+
+4. You can "drown" the `_FillValue` (i.e., set masked cells to near values) by:
+  - Using my `nn_interp.py` function: it's quite fast, and will fill a masked array with a the nearest neighbor *in terms of discrete coordinates* (i.e., array indices, *not* physical distance).
+  - Using the undocument `setmisstonn` operator in CDO: `cdo setmisstonn -setgrid,<cdo_griddes> masked.nc drowned`. In terms of code, it's straightforward and reliable, but this actually computes loads of distances (orthodromy) to get the physically nearest neighbor. So it's slow. But reliable. Just do it overnight!
+  - If your grid is regular (same lat/lon at each x/y), then better use `fillmiss` which fill masked values with linear reinterpolation.
+**Warning**: these drowning operations only work for array slices which have at least one unmasked cell in the (x,y) plane. If, on an array slice, all cells are masked, I think my `nn_interp.py` will just loop forever, and the CDO operators won't do anything. You probably need to do something else (like, replicate one vertical level to the other).
+
+5. If you want to perform the same interpolation a lot of times, it may be a good idea to generate a weight file, and then use it. The weight file is kind of the DNA of the interpolation. If CDO has it, it just uses the weight, and does a linear combination, which is much faster than interpolation without having the weights in advance. To do so, for bilinear interpolation:
+ ```
+ cdo genbil,<dest_cdo_griddes> -setgrid,<src_cdo_griddes> input.nc weights.nc
+ cdo remap,<dest_cdo_griddes>,weights.nc input.nc output.nc
+ ```
+ _Warning_: weight files depend on the mask of the input data in `input.nc`. If the input mask changes, e.g. if the data has a depth dimension along which the mask changes, this won't be useful.
+  
+6. Most times, conservative interpolation is "cleaner" but it's much harder to implement. In particular, you will need grid corners. Look at the CDO documentation and at the comments at the end of `get_cdo_griddes.sh`.
+
+7. Conservative interpolation from a masked input field is doable by:
+  - Putting zero's instead of _FillValue and removing the _FillValue attribute (this is an exception to point 3. above):
+     ```
+	 ncatted -a _FillValue,<some_var>,o,f,0. input.nc # store zero's on masked cells, but the netCDF will still consider them masked
+	 ncatted -a _FillValue,<some_var>,d,, input.nc # delete _FillValue attribute: the previous _FillValue will be there, hence, zero's
+	 ```
+   - Preparing a `float` or `double` precision input `mask` field, which is `0.` on masked cells and `1.` on unmasked cells.
+   - Perform conservative interpolation of the input data and of the mask.
+   - The result, on the destination grid and taking into account the mask, is `data/mask`. If `mask == 0`, this will look bad, but it's normal: you're looking at a destination masked cell (i.e. which has no overlap with any unmasked source cell).
+
+# CDO vertical interpolation
+
+Just the same as horizontal interpolation, CDO uses description files for describing a vertical axis. They can be generated from netCDF files like this:
+```
+cdo zaxisdes input.nc > description.zaxis
+```
+I never quite cornered what `input.nc` must exactly look like, but it should contain a "depth" field whose standard name is something like "depth", with informations on how it's stored. I know that raw 3D NEMO output fields work OK with this type of depth description:
+```
+	float deptht(deptht) ;
+		deptht:long_name = "Vertical T levels" ;
+		deptht:units = "m" ;
+		deptht:axis = "Z" ;
+		deptht:positive = "down" ;
+		deptht:bounds = "deptht_bounds" ;
+	float deptht_bounds(deptht, axis_nbounds) ;
+```