From 14ffd7d324496a66cc3d88dfd8b2b1b5ff765ba6 Mon Sep 17 00:00:00 2001
From: thopri <thopri@noc.ac.uk>
Date: Thu, 2 Apr 2020 14:52:59 +0100
Subject: [PATCH] updated fes tide code
---
inputs/namelist_remote.bdy | 9 ++--
pynemo/profile.py | 10 ++--
pynemo/tide/fes_extract_HC.py | 98 +++++++++++++++++++++++------------
pynemo/tide/nemo_bdy_tide3.py | 52 +++++++++++++------
4 files changed, 113 insertions(+), 56 deletions(-)
diff --git a/inputs/namelist_remote.bdy b/inputs/namelist_remote.bdy
index 9eeb1b6..97a61de 100644
--- a/inputs/namelist_remote.bdy
+++ b/inputs/namelist_remote.bdy
@@ -83,9 +83,12 @@
nn_month_end = 11 ! month end (default = 12 is years>1)
sn_dst_calendar = 'gregorian' ! output calendar format
nn_base_year = 1960 ! base year for time counter
- sn_tide_grid = './src_data/tide/grid_tpxo7.2.nc'
- sn_tide_h = './src_data/tide/h_tpxo7.2.nc'
- sn_tide_u = './src_data/tide/u_tpxo7.2.nc'
+ ! TPXO file locations
+ sn_tide_grid = '/Users/thopri/Projects/PyNEMO/DATA/TPXO/grid_tpxo7.2.nc'
+ sn_tide_h = '/Users/thopri/Projects/PyNEMO/DATA/TPXO/h_tpxo7.2.nc'
+ sn_tide_u = './Users/thopri/Projects/PyNEMO/DATA/TPXO/u_tpxo7.2.nc'
+ ! location of FES data
+ sn_tide_fes = '/Users/thopri/Projects/PyNEMO/DATA/FES/'
!------------------------------------------------------------------------------
! Additional parameters
diff --git a/pynemo/profile.py b/pynemo/profile.py
index a172225..83061b0 100644
--- a/pynemo/profile.py
+++ b/pynemo/profile.py
@@ -438,13 +438,15 @@ def process_bdy(setup_filepath=0, mask_gui=False):
if settings['tide']:
if settings['tide_model']=='tpxo':
- cons = tide.nemo_bdy_tpx7p2_rot(
+ cons = tide.nemo_bdy_tide_rot(
Setup, DstCoord, bdy_ind['t'], bdy_ind['u'], bdy_ind['v'],
- settings['clname'])
+ settings['clname'], settings['tide_model'])
elif settings['tide_model']=='fes':
- logger.error('Tidal model: %s, not yet implimented',
+ cons = tide.nemo_bdy_tide_rot(
+ Setup, DstCoord, bdy_ind['t'], bdy_ind['u'], bdy_ind['v'],
+ settings['clname'],settings['tide_model'])
+ logger.warning('Tidal model: %s, not yet properly implimented',
settings['tide_model'])
- return
else:
logger.error('Tidal model: %s, not recognised',
settings['tide_model'])
diff --git a/pynemo/tide/fes_extract_HC.py b/pynemo/tide/fes_extract_HC.py
index 7403781..747e043 100644
--- a/pynemo/tide/fes_extract_HC.py
+++ b/pynemo/tide/fes_extract_HC.py
@@ -12,8 +12,8 @@ from netCDF4 import Dataset
from scipy import interpolate
import numpy as np
-class FESExtract(object):
- """ This is TPXO model extract_hc.c implementation in python"""
+class FesExtract(object):
+ """ This is FES model extract_hc.c implementation in python adapted from tpxo_extract_HC.py"""
def __init__(self, settings, lat, lon, grid_type):
"""initialises the Extract of tide information from the netcdf
Tidal files"""
@@ -21,37 +21,69 @@ class FESExtract(object):
tide_model = 'fes'
if tide_model == 'fes': # Define stuff to generalise Tide model
- hRe_name = 'hRe'
- hIm_name = 'hIm'
- lon_z_name = 'lon_z'
- lat_z_name = 'lat_z'
- URe_name = 'URe'
- UIm_name = 'UIm'
- lon_u_name = 'lon_u'
- lat_u_name = 'lat_u'
- VRe_name = 'VRe'
- VIm_name = 'VIm'
- lon_v_name = 'lon_v'
- lat_v_name = 'lat_v'
- mz_name = 'mz'
- mu_name = 'mu'
- mv_name = 'mv'
- self.grid = Dataset(settings['tide_grid'])#../data/tide/grid_tpxo7.2.nc')
- #read the height_dataset file
- self.height_dataset = Dataset(settings['tide_h'])#../data/tide/h_tpxo7.2.nc')
- #read the velocity_dataset file
- self.velocity_dataset = Dataset(settings['tide_u'])#../data/tide/u_tpxo7.2.nc')
-
- height_z = self.grid.variables['hz']
- mask_z = self.grid.variables['mz']
- lon_z = self.height_dataset.variables[lon_z_name][:, 0]
- lat_z = self.height_dataset.variables[lat_z_name][0, :]
- lon_resolution = lon_z[1] - lon_z[0]
- data_in_km = 0 # added to maintain the reference to matlab tmd code
- # Pull out the constituents that are avaibable
- self.cons = []
- for ncon in range(self.height_dataset.variables['con'].shape[0]):
- self.cons.append(self.height_dataset.variables['con'][ncon, :].tostring().strip())
+
+ hRe_name = 'hRe'
+ hIm_name = 'hIm'
+ lon_z_name = 'lon'
+ lat_z_name = 'lat'
+ URe_name = 'URe'
+ UIm_name = 'UIm'
+ lon_u_name = 'lon'
+ lat_u_name = 'lat'
+ VRe_name = 'VRe'
+ VIm_name = 'VIm'
+ lon_v_name = 'lon'
+ lat_v_name = 'lat'
+ mz_name = 'mz'
+ mu_name = 'mu'
+ mv_name = 'mv'
+
+ #constituents = ['2N2', 'EPS2', 'J1', 'K1', 'K2', 'L2', 'LA2', 'M2', 'M3', 'M4', 'M6', 'M8', 'MF', 'MKS2',
+ #'MM', 'MN4', 'MS4', 'MSF', 'MSQM', 'MTM', 'MU2', 'N2', 'N4', 'NU2', 'O1', 'P1', 'Q1', 'R2',
+ #'S1', 'S2', 'S4', 'SA', 'SSA', 'T2']
+ constituents = ['M2','S2']
+
+ # extract lon and lat z data
+ lon_z = Dataset[settings['tide_fes']+constituents[0]+'_Z.nc'].variables['lon']
+ lat_z = Dataset[settings['tide_fes']+constituents[0]+'_Z.nc'].variabels['lat']
+ lon_resolution = lon_z[1] - lon_z[0]
+ data_in_km = 0 # added to maintain the reference to matlab tmd code
+ self.cons = constituents
+
+ #read and convert the height_dataset file to complex
+ for ncon in constituents:
+ amp = Dataset(settings['tide_fes']+constituents[ncon]+'_Z.nc').variables['amplitude'][:]
+ phase = Dataset(settings['tide_fes']+constituents[ncon]+'_Z.nc').variables['phase'][:]
+ hRe = amp*np.sin(phase)
+ hIm = amp*np.cos(phase)
+ self.height_dataset = {'hRe':hRe,'hIm':hIm}
+
+ #read and convert the velocity_dataset files to complex
+ for ncon in constituents:
+ amp = Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ua'][:]
+ phase = Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ug'][:]
+ URe = amp*np.sin(phase)
+ UIm = amp*np.cos(phase)
+ self.Uvelocity_dataset = {'URe':URe,'UIm':UIm}
+
+ for ncon in constituents:
+ amp = Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Va'][:]
+ phase = Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Vg'][:]
+ VRe = amp*np.sin(phase)
+ VIm = amp*np.cos(phase)
+ self.Vvelocity_dataset = {'VRe':VRe,'VIm':VIm}
+
+ # extract example amplitude grid and change NaNs to 0 (for land) and other values to 1 (for water)
+ mask_z = Dataset(settings['tide_fes']+constituents[ncon]+'_Z.nc').variables['amplitude'][:]
+ where_are_NaNs = np.isnan(mask_z)
+ where_no_NaNs = np.invert(np.isnan(mask_z))
+ mask_z[where_no_NaNs] = 1
+ mask_z[where_are_NaNs] = 0
+
+ # need to sort for FES
+ self.grid = Dataset(settings['tide_grid'])#../data/tide/grid_tpxo7.2.nc')
+ height_z = self.grid.variables['hz']
+
else:
print('Don''t know that tide model')
diff --git a/pynemo/tide/nemo_bdy_tide3.py b/pynemo/tide/nemo_bdy_tide3.py
index 7c9a148..4aa9222 100644
--- a/pynemo/tide/nemo_bdy_tide3.py
+++ b/pynemo/tide/nemo_bdy_tide3.py
@@ -8,6 +8,7 @@ Module to extract constituents for the input grid mapped onto output grid
# pylint: disable=no-name-in-module
import copy
from . import tpxo_extract_HC
+from . import fes_extract_HC
import numpy as np
from netCDF4 import Dataset
from pynemo import nemo_bdy_grid_angle
@@ -16,8 +17,8 @@ from pynemo.reader.factory import GetFile
import logging
-def nemo_bdy_tpx7p2_rot(setup, DstCoord, Grid_T, Grid_U, Grid_V, comp):
- """ TPXO Global Tidal model interpolation including rotation grid"""
+def nemo_bdy_tide_rot(setup, DstCoord, Grid_T, Grid_U, Grid_V, comp,tide_model):
+ """ Global Tidal model interpolation including rotation grid"""
key_transport = 0 # compute the velocities from transport
numharm = len(comp)
logger = logging.getLogger(__name__)
@@ -52,16 +53,25 @@ def nemo_bdy_tpx7p2_rot(setup, DstCoord, Grid_T, Grid_U, Grid_V, comp):
#convert the U-longitudes into the TMD conventions (0/360E)
dst_lon[dst_lon < 0.0] = dst_lon[dst_lon < 0.0]+360.0
+ if tide_model == 'tpxo':
+ tpxo_ux = tpxo_extract_HC.TpxoExtract(setup.settings, dst_lat, dst_lon, Grid_U.grid_type)
+ tpxo_vx = tpxo_extract_HC.TpxoExtract(setup.settings, dst_lat, dst_lon, Grid_V.grid_type)
- tpxo_ux = tpxo_extract_HC.TpxoExtract(setup.settings, dst_lat, dst_lon, Grid_U.grid_type)
- tpxo_vx = tpxo_extract_HC.TpxoExtract(setup.settings, dst_lat, dst_lon, Grid_V.grid_type)
+ ampuX = tpxo_ux.amp
+ phauX = tpxo_ux.gph
+ ampvX = tpxo_vx.amp
+ phavX = tpxo_vx.gph
- ampuX = tpxo_ux.amp
- phauX = tpxo_ux.gph
- ampvX = tpxo_vx.amp
- phavX = tpxo_vx.gph
+ if tide_model == 'fes':
+ fes_uy = fes_extract_HC.FesExtract(setup.settings, dst_lat, dst_lon, Grid_U.grid_type)
+ fes_vy = fes_extract_HC.FesExtract(setup.settings, dst_lat, dst_lon, Grid_V.grid_type)
- #check if ux data are missing
+ ampuY = fes_uy.amp
+ phauY = fes_uy.gph
+ ampvY = fes_vy.amp
+ phavY = fes_vy.gph
+
+ #check if ux data are missing
ind = np.where((np.isnan(ampuX)) | (np.isnan(phauX)))
if ind[0].size > 0:
logger.warning('Missing zonal velocity along the x open boundary')
@@ -83,13 +93,23 @@ def nemo_bdy_tpx7p2_rot(setup, DstCoord, Grid_T, Grid_U, Grid_V, comp):
#convert the U-longitudes into the TMD conventions (0/360E)
dst_lon[dst_lon < 0.0] = dst_lon[dst_lon < 0.0]+360.0
- tpxo_uy = tpxo_extract_HC.TpxoExtract(setup.settings, dst_lat, dst_lon, Grid_U.grid_type)
- tpxo_vy = tpxo_extract_HC.TpxoExtract(setup.settings, dst_lat, dst_lon, Grid_V.grid_type)
-
- ampuY = tpxo_uy.amp
- phauY = tpxo_uy.gph
- ampvY = tpxo_vy.amp
- phavY = tpxo_vy.gph
+ if tide_model == 'tpxo':
+ tpxo_uy = tpxo_extract_HC.TpxoExtract(setup.settings, dst_lat, dst_lon, Grid_U.grid_type)
+ tpxo_vy = tpxo_extract_HC.TpxoExtract(setup.settings, dst_lat, dst_lon, Grid_V.grid_type)
+
+ ampuY = tpxo_uy.amp
+ phauY = tpxo_uy.gph
+ ampvY = tpxo_vy.amp
+ phavY = tpxo_vy.gph
+
+ if tide_model == 'fes':
+ fes_uy = fes_extract_HC.FesExtract(setup.settings, dst_lat, dst_lon, Grid_U.grid_type)
+ fes_vy = fes_extract_HC.FesExtract(setup.settings, dst_lat, dst_lon, Grid_V.grid_type)
+
+ ampuY = fes_uy.amp
+ phauY = fes_uy.gph
+ ampvY = fes_vy.amp
+ phavY = fes_vy.gph
#check if ux data are missing
ind = np.where((np.isnan(ampuY)) | (np.isnan(phauY)))
--
GitLab