Skip to content

Commit

Permalink
Read unstructured domain resample (#18)
Browse files Browse the repository at this point in the history
* Initial commit of unstructured mesh regrided to rectilinear grid

* Added more comprehensive command line parsing to get_netcdf.py

* Add default values to help text

* Work out mesh topology variable name if not given. Improve help text

* Fix using global var instead of local var bug. Put main code into function so it can be imported and called

* Add code to make tests work

* Changes requested in pull request, delete unused python script and tidy up shared_testing.py

* Add extra comments/docstrings

* Remove XIOS3 skip test

* Increase relative tolerance value so there are no known failures

* Replace Dimension instances with Dimension names in createVariable as only versions >=1.6.0 of python netCDF4 support using Dimension instances

* Add read_access=.true. to xml, seems to be needed by XIOS3 but not XIOS2

* Replace mesh_C12.nc with cdl version + code changes to use cdl file

* Add back XIOS3 skip test
  • Loading branch information
jeff-cole authored May 22, 2024
1 parent 803eeb2 commit 10cde6a
Show file tree
Hide file tree
Showing 13 changed files with 6,013 additions and 5 deletions.
124 changes: 124 additions & 0 deletions xios_examples/dataFunc.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,124 @@
import numpy as np

class dataFunc:
"""
Class containing functions which generate analytical data values for an input of
latitude and longitude arrays.
Functions sinusiod,harmonic,vortex,gulfstream are taken from the paper
`Benchmarking Regridding Libraries Used in Earth System Modelling`, see
https://www.mdpi.com/2297-8747/27/2/31
"""

def func_sinusiod(self, latarr, lonarr):

length = 1.2*np.pi
conv = np.pi/180.0
coef = 2.0
coefmult = 1.0

data = np.array(coefmult*(coef - np.cos( np.pi*(np.arccos( np.cos(lonarr*conv)*np.cos(latarr*conv) )/length))), dtype=np.float64)

return data

def func_harmonic(self, latarr, lonarr):

conv = np.pi/180.0

data = np.array(2.0 + (np.sin(2.0*latarr*conv)**16)*np.cos(16.0*lonarr*conv), dtype=np.float64)

return data

def func_vortex(self, latarr, lonarr):

lon0 = 5.5
lat0 = 0.2
r0 = 3.0
d = 5.0
t = 6.0
conv = np.pi/180.0
sinc = np.sin(lat0)
cosc = np.cos(lat0)

# Find the rotated Longitude and Latitude of a point on a sphere
# with pole at (lon0, lat0)
cost = np.cos(latarr*conv)
sint = np.sin(latarr*conv)

trm = cost * np.cos(lonarr*conv-lon0)
x = sinc * trm - cosc*sint
y = cost * np.sin(lonarr*conv-lon0)
z = sinc * sint + cosc*trm

lon = np.arctan2(y, x)
lon = np.where(lon < 0.0, lon+2.0*np.pi, lon)
lat = np.arcsin(z)

rho = r0 * np.cos(lat)
vt = 3.0 * np.sqrt(3.0)/2.0/np.cosh(rho)/np.cosh(rho)*np.tanh(rho)
omega = np.where(rho == 0.0, 0.0, vt/rho)

data = np.array(2.0*(1.0+np.tanh(rho/d * np.sin(lon-omega*t))), dtype=np.float64)

return data

def func_gulfstream(self, latarr, lonarr):

# Parameters for analytical function
coef = 2.0
length = 1.2*np.pi
conv = np.pi/180.0
gf_coef = 1.0 # Coefficient for Gulf Stream term (0.0 = no Gulf Stream)
gf_ori_lon = -80.0 # Origin of the Gulf Stream (longitude in deg)
gf_ori_lat = 25.0 # Origin of the Gulf Stream (latitude in deg)
gf_end_lon = -1.8 # End point of the Gulf Stream (longitude in deg)
gf_end_lat = 50.0 # End point of the Gulf Stream (latitude in deg)
gf_dmp_lon = -25.5 # Point of the Gulf Stream decrease (longitude in deg)
gf_dmp_lat = 55.5 # Point of the Gulf Stream decrease (latitude in deg)

dr0 = np.sqrt(((gf_end_lon - gf_ori_lon)*conv)**2 +
((gf_end_lat - gf_ori_lat)*conv)**2)

dr1 = np.sqrt(((gf_dmp_lon - gf_ori_lon)*conv)**2 +
((gf_dmp_lat - gf_ori_lat)*conv)**2)

# Original OASIS fcos analytical test function
fnc_ana = (coef-np.cos(np.pi*(np.arccos(np.cos(latarr*conv)*np.cos(lonarr*conv))/length)))
gf_per_lon = lonarr
gf_per_lon = np.where(gf_per_lon > 180.0, gf_per_lon-360.0, gf_per_lon)
gf_per_lon = np.where(gf_per_lon < -180.0, gf_per_lon+360.0, gf_per_lon)
dx = (gf_per_lon - gf_ori_lon)*conv
dy = (latarr - gf_ori_lat)*conv
dr = np.sqrt(dx*dx + dy*dy)
dth = np.arctan2(dy, dx)
dc = 1.3*gf_coef
dc = np.where(dr > dr0, 0.0, dc)
dc = np.where(dr > dr1, dc * np.cos(np.pi*0.5*(dr-dr1)/(dr0-dr1)), dc)
data = np.array(fnc_ana + (np.maximum(1000.0*np.sin(0.4*(0.5*dr+dth)+0.007*np.cos(50.0*dth) +
0.37*np.pi),999.0)-999.0)*dc, dtype=np.float64)

return data

def func_cossin(self, latarr, lonarr):

length = 1.0*np.pi
conv = np.pi/180.
coef = 21.0
coefmult = 3.846 * 20.0

data = np.array(coefmult*(coef - np.cos( np.pi*(np.arccos( np.cos(lonarr*conv)*np.cos(latarr*conv) )/length)) *
np.sin( np.pi*(np.arcsin( np.sin(lonarr*conv)*np.sin(latarr*conv) )/length))), dtype=np.float64)

return data

def get_funclist(self):

funclist = [func.removeprefix('func_') for func in dir(self) if callable(getattr(self, func)) and func.startswith('func_')]

return funclist

def get_func(self, name: str):

do = f"func_{name}"
if hasattr(self, do) and callable(func := getattr(self, do)):
return func
Loading

0 comments on commit 10cde6a

Please sign in to comment.