Commit 03b45075a0d by Pierre Augier

Start to add some functions.

parent 008ef33b14e
......@@ -7,6 +7,68 @@
"""
import os
import numpy as np
class DataSetGCM(object):
"""Represent a dataset."""
def __init__(self, name=None):
self.name = name
self.path_result = ''
self.path_data = ''
self.plevs = None
def print_info(self):
print('\nInformation: \"' + self.name + '\"'
'nlon = {0}, nlat = {1}'.format(self.nlon, self.nlat) +
'1 point every {0:6.2f} km'.format(
2*np.pi*r_a_real/self.nlon/1000))
if hasattr(self, 'nz'):
print('{0} levels in the vertical'.format(self.nz))
print('')
if hasattr(self, 'esh'):
print('Info library shtns (spherical harmonic transforms):')
self.esh.print_info()
print('')
def load_spectral3dvar(self, key, instant):
raise NotImplementedError
def load_spatial3dvar(self, key, instant):
raise NotImplementedError
def load_spectral2dvar(self, key, instant):
raise NotImplementedError
def load_spatial2dvar(self, key, instant):
raise NotImplementedError
def compute_timeaverages_beta(self):
raise NotImplementedError
def compute_timeaverages_gamma(self):
raise NotImplementedError
def compute_hmean(self, f_xy):
"""Compute the horizontal mean."""
delta = 360./self.nb_lon
return (delta*np.pi/180)**2 * np.sum(f_xy*self.cosLATS) / (4*np.pi)
def compute_hmean_representative(self, field, ip=None, beta=None):
"""Compute the horizontal "representative" mean."""
if beta is None:
if 'beta3d' not in dir(self):
raise ValueError(
'first run sim.load_beta_ps_mean() ???')
if ip < self.nb_lev_beta:
result = self.mean_field(
field*self.beta3d[ip])/self.beta_mean1d[ip]
elif ip >= self.nb_lev_beta:
result = self.mean_field(field)
else:
result = self.mean_field(field*beta)/self.mean_field(beta)
return result
"""Define some physical constants
=================================
"""
c_p = 1006. # (J/kg/K)
R = 287. # (J/kg/K)
g = 9.80665 # (m/s^2)
Gamma_dah = g/c_p
khi = R/c_p
......@@ -5,25 +5,29 @@
:members:
:private-members:
.. autoclass:: Earth
:members:
:private-members:
.. autoclass:: Mars
:members:
:private-members:
"""
from math import pi
class Planet(object):
"""Represent a planet."""
"""Represent a planet.
...
"""
class Earth(Planet):
"""Represent the Earth."""
def __init__(self, name, radius, Omega, pressure_m):
self.name = 'Earth'
self.radius = radius
self.Omega = Omega
self.pressure_m = pressure_m
class Mars(Planet):
"""Represent Mars."""
earth = Planet(
'Earth',
radius=6367470., # Earth radius (meters)
Omega=2*pi/86164, # Rotation frequency (Hz)
pressure_m=1013. # Mean pressure at sea level (hPa)
)
......@@ -121,18 +121,18 @@ class EasySHT(object):
- lat denotes latitude
- lon denotes longitude
- kh denotes horizontal wavenumber
- uu denotes longitudinal velocity
- vv denotes meridional velocity (the sign depends on the used base)
- u denotes longitudinal velocity
- v denotes meridional velocity (the sign depends on the used base)
- hdiv denotes horizontal divergence
- hrot denotes vertical vorticity (curl on the sphere)
- grad denotes the horizontal gradient
Variables and functions about spectral space array:
lmax, mmax, mres, nlm, idx_lm(), l_idx, m_idx, l2_idx
lmax, mmax, mres, nb_lm, idx_lm(), l_idx, m_idx, l2_idx
Variables about grid and real space array:
nlat, nlon, lats, lons, sin_lats, LATS, LONS
nb_lat, nb_lon, lats, lons, sin_lats, LATS, LONS
Variables for spectra:
l2_l, kh_l
......@@ -192,7 +192,7 @@ class EasySHT(object):
self.l_idx = self.sh.l
self.l2_idx = self.l_idx*(self.l_idx+1)
self.m_idx = self.sh.m
self.nlm = self.sh.nlm
self.nb_lm = self.sh.nlm
self.delta = 360./(self.nlon*self.mres)
# create arrays 2D lats et lons
......@@ -203,7 +203,7 @@ class EasySHT(object):
self.l2_l = lrange*(lrange+1)
self.kh_l = np.sqrt(self.l2_l)/self.radius
self.complex64_save_netCFD = np.dtype(
self._complex64_save_netCFD = np.dtype(
[('real', np.float32), ('imag', np.float32)])
def idx_lm(self, l, m):
......@@ -218,13 +218,13 @@ class EasySHT(object):
def create_array_SH(self, value=None, dtype=complex):
"""Create an array representing a field in spectral space."""
if value is None:
field_lm = np.empty(self.nlm, dtype)
field_lm = np.empty(self.nb_lm, dtype)
elif value=='rand':
field_lm = np.random.randn(self.nlm)+1.j*np.random.randn(self.nlm)
field_lm = np.random.randn(self.nb_lm)+1.j*np.random.randn(self.nb_lm)
elif value==0:
field_lm = np.zeros(self.nlm, dtype)
field_lm = np.zeros(self.nb_lm, dtype)
else:
field_lm = value*np.ones(self.nlm, dtype)
field_lm = value*np.ones(self.nb_lm, dtype)
return field_lm
def create_array_spat(self, value=None):
......@@ -241,7 +241,7 @@ class EasySHT(object):
return field
def convert2complex64_save_netCFD(self, f_lm):
result = np.empty(f_lm.shape, self.complex64_save_netCFD)
result = np.empty(f_lm.shape, self._complex64_save_netCFD)
result['real'] = f_lm.real
result['imag'] = f_lm.imag
return result
......@@ -309,10 +309,10 @@ class EasySHT(object):
+''.format((t2-t1)/nb_sht))
# functions for 2D vectorial spherical harmonic transforms
def uuvv_from_hdivrotSH(self, hdiv_lm, hrot_lm, uu=None, vv=None,
PRINT_TIME=False):
def uv_from_hdivrotSH(self, hdiv_lm, hrot_lm, uu=None, vv=None,
PRINT_TIME=False):
"""
uu vv from h, div, rot (uu and vv are overwritten)
u, v from h, div, rot (u and v are overwritten)
"""
if PRINT_TIME:
t1 = time()
......@@ -333,10 +333,10 @@ class EasySHT(object):
'{:4.3f} s'.format(t2-t1))
return uu, vv
def uuvv_from_uDuRSH(self, uD_lm, uR_lm, uu=None, vv=None,
PRINT_TIME=False):
def uv_from_uDuRSH(self, uD_lm, uR_lm, uu=None, vv=None,
PRINT_TIME=False):
"""
uu vv from uD, uR (uu and vv are overwritten)
u, v from uD, uR (uu and vv are overwritten)
"""
if PRINT_TIME:
t1 = time()
......@@ -356,8 +356,8 @@ class EasySHT(object):
def hdivrotSH_from_uuvv(self, uu, vv, hdiv_lm=None, hrot_lm=None,
PRINT_TIME=False):
def hdivrotSH_from_uv(self, uu, vv, hdiv_lm=None, hrot_lm=None,
PRINT_TIME=False):
"""
hdivrotSH_from_uuvv
(div_lm and rot_lm are overwritten)
......@@ -409,7 +409,7 @@ class EasySHT(object):
t1 = time()
if gradf_lon is None:
gradf_lon = self.create_array_spat(0)
gradf_lat = self.create_array_spat(0) # becareful bug if not 0!!!
gradf_lat = self.create_array_spat(0) # becareful bug if not 0!!!
# We do not use SHsph_to_spat() because it seems that there is a problem
#### self.sh.SHsph_to_spat(f_lm, gradf_lat, gradf_lon)
# instead we use SHsphtor_to_spat(...) with tor_lm= zeros_lm
......@@ -455,7 +455,7 @@ class EasySHT(object):
t1 = time()
nvert = f3D.shape[0]
if f_lm3D is None:
f_lm3D = np.empty([nvert, self.nlm], dtype)
f_lm3D = np.empty([nvert, self.nb_lm], dtype)
for iz in range(nvert):
f_lm3D[iz] = self.SH_from_spat(f3D[iz])
if PRINT_TIME:
......@@ -491,15 +491,15 @@ class EasySHT(object):
def _spectrum_from_array_deSH(self, array_deSH):
""" compute spectrum(l) from array_deSH(ilm)"""
spectrum = np.zeros(self.lmax+1)
for ilm in range(0, self.nlm):
for ilm in range(0, self.nb_lm):
spectrum[self.l_idx[ilm]] += array_deSH[ilm]
return spectrum
def _array_deSH_from_SH(self, field_lm, key_field):
""" compute the array_deSH (density of energy) from an field_lm"""
if key_field[:2]=='uu' or key_field[:2]=='vv':
if key_field[:1]=='u' or key_field[:1]=='v':
array_deSH = abs(field_lm)**2/2.
elif key_field[:2]=='TT' or key_field[:2]=='ps' or key_field[:2]=='oo':
elif key_field[:1]=='T' or key_field[:2]=='ps' or key_field[:1]=='o':
array_deSH = abs(field_lm)**2/2.
elif key_field[:4]=='beta':
array_deSH = self.create_array_SH(0., float)
......@@ -508,7 +508,7 @@ class EasySHT(object):
array_deSH = self.l2_idx* (abs(field_lm)**2)/2
elif key_field[:4]=='hdiv' or key_field[:4]=='hrot':
array_deSH = self.create_array_SH(0., float)
COND = self.l2_idx>0
COND = self.l2_idx > 0
array_deSH[COND] = self.radius**2/self.l2_idx[COND]*abs(
field_lm[COND])**2/2
else:
......@@ -522,29 +522,32 @@ class EasySHT(object):
spectrum = self._spectrum_from_array_deSH(array_deSH)
return spectrum
def cospectrum_from_2fieldsSH(self, ff_lm, gg_lm):
""" compute cospectrum(l) from ff_lm(ilm) and gg_lm(ilm)"""
def cospectrum_from_2fieldsSH(self, f_lm, g_lm):
""" compute cospectrum(l) from f_lm(ilm) and g_lm(ilm)"""
cospectrum = np.zeros(self.lmax+1)
array_deSH = ff_lm.conjugate()*gg_lm + ff_lm*gg_lm.conjugate()
array_deSH = f_lm.conjugate() * g_lm + f_lm * g_lm.conjugate()
array_deSH = array_deSH.real
array_deSH[self.m_idx==0] = array_deSH[self.m_idx==0]/2
for ilm in range(0, self.nlm):
for ilm in range(self.nb_lm):
cospectrum[self.l_idx[ilm]] += array_deSH[ilm]
return cospectrum
def cospectrum_from_2vectorsSH(self, ff_lon_lm, ff_lat_lm,
gg_lon_lm, gg_lat_lm):
def cospectrum_from_2vectorsSH(self, f_lon_lm, f_lat_lm,
g_lon_lm, g_lat_lm):
""" compute cospectrum(l)..."""
cospectrum = np.zeros(self.lmax+1)
array_deSH = (
ff_lon_lm.conjugate()*gg_lon_lm + ff_lon_lm*gg_lon_lm.conjugate() +
ff_lat_lm.conjugate()*gg_lat_lm + ff_lat_lm*gg_lat_lm.conjugate())
f_lon_lm.conjugate() * g_lon_lm + f_lon_lm * g_lon_lm.conjugate() +
f_lat_lm.conjugate() * g_lat_lm + f_lat_lm * g_lat_lm.conjugate())
array_deSH = array_deSH.real
array_deSH[self.m_idx==0] = array_deSH[self.m_idx==0]/2
for ilm in range(self.nlm):
for ilm in range(self.nb_lm):
cospectrum[self.l_idx[ilm]] += array_deSH[ilm]
return cospectrum
def cospectrum_from_2divrotSH(self):
raise NotImplementedError
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment