Commit c81a69dd authored by Dion Häfner's avatar Dion Häfner

pull in changes from mobi branch

parent e73ca620
...@@ -184,3 +184,85 @@ def tempsalt_sources(vs): ...@@ -184,3 +184,85 @@ def tempsalt_sources(vs):
vs.int_drhodS[..., vs.tau] * vs.salt_source) vs.int_drhodS[..., vs.tau] * vs.salt_source)
vs.P_diss_sources[...] = 0. vs.P_diss_sources[...] = 0.
dissipation_on_wgrid(vs, vs.P_diss_sources, aloc=aloc) dissipation_on_wgrid(vs, vs.P_diss_sources, aloc=aloc)
@veros_method
def biharmonic(vs, tr, f, dtr):
"""
Biharmonic mixing of tracer tr, results saved as dtr
This is essentially just a copy of tempsalt_biharmonic generalized
"""
flux_east = np.zeros_like(tr, dtype=vs.default_float_type)
flux_north = np.zeros_like(tr, dtype=vs.default_float_type)
del2 = np.zeros((vs.nx + 4, vs.ny + 4, vs.nz), dtype=vs.default_float_type)
flux_east[:-1, :, :] = -f * (tr[1:, :, :] - tr[:-1, :, :]) \
/ (vs.cost[np.newaxis, :, np.newaxis] * vs.dxu[:-1, np.newaxis, np.newaxis]) \
* vs.maskU[:-1, :, :]
flux_north[:, :-1, :] = -f * (tr[:, 1:, :] - tr[:, :-1, :]) \
/ vs.dyu[np.newaxis, :-1, np.newaxis] * vs.maskV[:, :-1, :] \
* vs.cosu[np.newaxis, :-1, np.newaxis]
flux_east[:, -1, :] = 0.
flux_north[:, -1, :] = 0.
del2[1:, 1:, :] = vs.maskT[1:, 1:, :] * (flux_east[1:, 1:, :] - flux_east[:-1, 1:, :]) \
/ (vs.cost[np.newaxis, 1:, np.newaxis] * vs.dxt[1:, np.newaxis, np.newaxis]) \
+ (flux_north[1:, 1:, :] - flux_north[1:, :-1, :]) \
/ (vs.cost[np.newaxis, 1:, np.newaxis] * vs.dyt[np.newaxis, 1:, np.newaxis])
utilities.enforce_boundaries(vs, del2)
flux_east[:-1, :, :] = f * (del2[1:, :, :] - del2[:-1, :, :]) \
/ (vs.cost[np.newaxis, :, np.newaxis] * vs.dxu[:-1, np.newaxis, np.newaxis]) \
* vs.maskU[:-1, :, :]
flux_north[:, :-1, :] = f * (del2[:, 1:, :] - del2[:, :-1, :]) \
/ vs.dyu[np.newaxis, :-1, np.newaxis] * vs.maskV[:, :-1, :] \
* vs.cosu[np.newaxis, :-1, np.newaxis]
flux_east[-1, :, :] = 0.
flux_north[:, -1, :] = 0.
dtr[1:, 1:, :] = vs.maskT[1:, 1:, :] * (flux_east[1:, 1:, :] - flux_east[:-1, 1:, :]) \
/ (vs.cost[np.newaxis, 1:, np.newaxis] * vs.dxt[1:, np.newaxis, np.newaxis]) \
+ (flux_north[1:, 1:, :] - flux_north[1:, :-1, :]) \
/ (vs.cost[np.newaxis, 1:, np.newaxis] * vs.dyt[np.newaxis, 1:, np.newaxis])
dtr[...] *= vs.maskT
# if vs.enable_conserve_energy: TODO should we do something here?
@veros_method
def horizontal_diffusion(vs, tr, dtr_hmix):
"""
Diffusion of tracer tr
"""
# aloc = np.zeros((vs.nx + 4, vs.ny + 4, vs.nz), dtype=vs.default_float_type)
flux_east = np.zeros_like(tr, dtype=vs.default_float_type)
flux_north = np.zeros_like(tr, dtype=vs.default_float_type)
# horizontal diffusion of tracer
flux_east[:-1, :, :] = vs.K_h * (tr[1:, :, :] - tr[:-1, :, :]) \
/ (vs.cost[np.newaxis, :, np.newaxis] * vs.dxu[:-1, np.newaxis, np.newaxis])\
* vs.maskU[:-1, :, :]
flux_east[-1, :, :] = 0.
flux_north[:, :-1, :] = vs.K_h * (tr[:, 1:, :] - tr[:, :-1, :]) \
/ vs.dyu[np.newaxis, :-1, np.newaxis] * vs.maskV[:, :-1, :]\
* vs.cosu[np.newaxis, :-1, np.newaxis]
flux_north[:, -1, :] = 0.
if vs.enable_hor_friction_cos_scaling:
flux_east[...] *= vs.cost[np.newaxis, :, np.newaxis] ** vs.hor_friction_cosPower
flux_north[...] *= vs.cosu[np.newaxis, :, np.newaxis] ** vs.hor_friction_cosPower
dtr_hmix[1:, 1:, :] = ((flux_east[1:, 1:, :] - flux_east[:-1, 1:, :])
/ (vs.cost[np.newaxis, 1:, np.newaxis] * vs.dxt[1:, np.newaxis, np.newaxis])
+ (flux_north[1:, 1:, :] - flux_north[1:, :-1, :])
/ (vs.cost[np.newaxis, 1:, np.newaxis] * vs.dyt[np.newaxis, 1:, np.newaxis]))\
* vs.maskT[1:, 1:, :]
# if vs.enable_conserve_energy:
# dissipation_on_wgrid(vs, vs.P_diss_hmix, int_drhodX=vs.int_drhodT[..., vs.tau])
...@@ -97,6 +97,45 @@ def _calc_implicit_part(vs, tr): ...@@ -97,6 +97,45 @@ def _calc_implicit_part(vs, tr):
tr[2:-2, 2:-2, :, vs.taup1] = utilities.where(vs, water_mask, sol, tr[2:-2, 2:-2, :, vs.taup1]) tr[2:-2, 2:-2, :, vs.taup1] = utilities.where(vs, water_mask, sol, tr[2:-2, 2:-2, :, vs.taup1])
@veros_method
def isoneutral_diffusion_decoupled(vs, tr, dtracer_iso, iso=True, skew=False):
"""
Like isoneutral_diffusion but for general tracers
"""
if iso:
K_iso = vs.K_iso
else:
K_iso = np.zeros_like(vs.K_iso)
if skew:
K_skew = vs.K_gm
else:
K_skew = np.zeros_like(vs.K_gm)
_calc_tracer_fluxes(vs, tr, K_iso, K_skew)
"""
add explicit part
"""
aloc = _calc_explicit_part(vs)
dtracer_iso[...] += aloc[...]
tr[2:-2, 2:-2, :, vs.taup1] += vs.dt_tracer * aloc[2:-2, 2:-2, :]
"""
add implicit part
"""
if iso:
aloc[...] = tr[:, :, :, vs.taup1]
_calc_implicit_part(vs, tr)
dtracer_iso[...] += (tr[:, :, :, vs.taup1] - aloc) / vs.dt_tracer
"""
dissipation by isopycnal mixing
"""
# TODO add dissipation by isopycnal mixing
# NOTE drhodS, drhodT are calculated in isoneutral.py. The rest in there is additional stuff
@veros_method @veros_method
def isoneutral_diffusion(vs, tr, istemp, iso=True, skew=False): def isoneutral_diffusion(vs, tr, istemp, iso=True, skew=False):
""" """
......
from veros.setup.bgc_global_4deg.bgc_global_four_degree import GlobalFourDegreeBGC
forcing:
url: https://sid.erda.dk/share_redirect/gsdZADr8to/global_4deg/forcing_4deg_global.nc
md5: ef3be0a58782771c8ee5a6d0206b87f5
ecmwf:
url: https://sid.erda.dk/share_redirect/gsdZADr8to/global_4deg/ecmwf_4deg_monthly_nc4.nc
md5: d1b4e0e199d7a5883cf7c88d3d6bcb27
corev2:
url: https://sid.erda.dk/share_redirect/gsdZADr8to/COREv2_NYF/CNYFv2/monthly_average/ncar_rad.15JUNE2009.monthly_avg_4deg.nc
md5: 328a1bcf37b2b14671c8526e462cc1c1
#!/usr/bin/env python
import os
import h5netcdf
import ruamel.yaml as yaml
from veros import VerosSetup, veros_method, distributed
from veros.variables import Variable
import veros.tools
BASE_PATH = os.path.dirname(os.path.realpath(__file__))
DATA_FILES = veros.tools.get_assets(
'global_4deg',
os.path.join(BASE_PATH, 'assets.yml')
)
class GlobalFourDegreeBGC(VerosSetup):
"""Global 4 degree model with 15 vertical levels and biogeochemistry.
"""
@veros_method
def set_parameter(self, vs):
vs.identifier = '4deg'
vs.nx, vs.ny, vs.nz = 90, 40, 15
vs.dt_mom = 1800.0
vs.dt_tracer = 86400.0
vs.dt_bio = vs.dt_tracer // 4
vs.runlen = 0.
vs.trcmin = 0
vs.enable_npzd = True
vs.enable_carbon = True
with open(os.path.join(BASE_PATH, 'npzd.yml')) as yaml_file:
cfg = yaml.safe_load(yaml_file)['npzd']
vs.npzd_selected_rules = cfg['selected_rules']
vs.remineralization_rate_detritus = 0.09 / 86400
vs.bbio = 1.038
vs.cbio = 1.0
vs.maximum_growth_rate_phyto = 0.23 / 86400
vs.maximum_grazing_rate = 0.13 / 86400
vs.fast_recycling_rate_phytoplankton = 0.025 / 86400
vs.specific_mortality_phytoplankton = 0.035 / 86400
vs.quadric_mortality_zooplankton = 0.06 / 86400
vs.zooplankton_growth_efficiency = 0.60
vs.assimilation_efficiency = 0.5
vs.wd0 = 2 / 86400
vs.mwz = 1000
vs.mw = 0.02 / 86400
vs.dcaco3 = 2500
vs.coord_degree = True
vs.enable_cyclic_x = True
vs.congr_epsilon = 1e-8
vs.congr_max_iterations = 20000
vs.enable_neutral_diffusion = True
vs.K_iso_0 = 1000.0
vs.K_iso_steep = 500.0
vs.iso_dslope = 0.001
vs.iso_slopec = 0.001
vs.enable_skew_diffusion = True
vs.enable_hor_friction = True
vs.A_h = (4 * vs.degtom)**3 * 2e-11
vs.enable_hor_friction_cos_scaling = True
vs.hor_friction_cosPower = 1
vs.enable_implicit_vert_friction = True
vs.enable_tke = True
vs.c_k = 0.1
vs.c_eps = 0.7
vs.alpha_tke = 30.0
vs.mxl_min = 1e-8
vs.tke_mxl_choice = 2
vs.kappaM_min = 2e-4
vs.kappaH_min = 7e-5 # usually 2e-5
vs.enable_Prandtl_tke = False
vs.enable_kappaH_profile = True
# eke
self.K_gm_0 = 1000.0
self.enable_eke = False
self.eke_k_max = 1e4
self.eke_c_k = 0.4
self.eke_c_eps = 0.5
self.eke_cross = 2.
self.eke_crhin = 1.0
self.eke_lmin = 100.0
self.enable_eke_superbee_advection = False
self.enable_eke_isopycnal_diffusion = False
# idemix
self.enable_idemix = False
self.enable_idemix_hor_diffusion = False
self.enable_eke_diss_surfbot = False
self.eke_diss_surfbot_frac = 0.2 # fraction which goes into bottom
self.enable_idemix_superbee_advection = False
vs.eq_of_state_type = 5
# custom variables
vs.nmonths = 12
vs.variables.update(
sss_clim=Variable('sss_clim', ('xt', 'yt', 'nmonths'), '', '', time_dependent=False),
sst_clim=Variable('sst_clim', ('xt', 'yt', 'nmonths'), '', '', time_dependent=False),
qnec=Variable('qnec', ('xt', 'yt', 'nmonths'), '', '', time_dependent=False),
qnet=Variable('qnet', ('xt', 'yt', 'nmonths'), '', '', time_dependent=False),
taux=Variable('taux', ('xt', 'yt', 'nmonths'), '', '', time_dependent=False),
tauy=Variable('tauy', ('xt', 'yt', 'nmonths'), '', '', time_dependent=False),
)
@veros_method
def _read_forcing(self, vs, var):
with h5netcdf.File(DATA_FILES['forcing'], 'r') as infile:
var_obj = infile.variables[var]
return np.array(var_obj, dtype=str(var_obj.dtype)).T
@veros_method
def set_grid(self, vs):
ddz = np.array([50., 70., 100., 140., 190., 240., 290., 340.,
390., 440., 490., 540., 590., 640., 690.])
vs.dzt[:] = ddz[::-1]
vs.dxt[:] = 4.0
vs.dyt[:] = 4.0
vs.y_origin = -78.0
vs.x_origin = 4.0
idx_global, _ = distributed.get_chunk_slices(vs, ('xt', 'yt', 'zt'))
with h5netcdf.File(DATA_FILES['corev2'], 'r') as infile:
vs.swr_initial = infile.variables['SWDN_MOD'][idx_global[::-1]].T
@veros_method
def set_coriolis(self, vs):
vs.coriolis_t[...] = 2 * vs.omega * np.sin(vs.yt[np.newaxis, :] / 180. * vs.pi)
@veros_method(dist_safe=False, local_variables=[
'kbot'
])
def set_topography(self, vs):
bathymetry_data = self._read_forcing(vs, 'bathymetry')
salt_data = self._read_forcing(vs, 'salinity')[:, :, ::-1]
mask_salt = salt_data == 0.
vs.kbot[2:-2, 2:-2] = 1 + np.sum(mask_salt.astype(np.int), axis=2)
mask_bathy = bathymetry_data == 0
vs.kbot[2:-2, 2:-2][mask_bathy] = 0
vs.kbot[vs.kbot == vs.nz] = 0
@veros_method(dist_safe=False, local_variables=[
'taux', 'tauy', 'qnec', 'qnet', 'sss_clim', 'sst_clim',
'temp', 'salt', 'taux', 'tauy', 'area_t', 'maskT',
'forc_iw_bottom', 'forc_iw_surface', 'zw', 'maskT',
'phytoplankton', 'zooplankton', 'detritus', 'po4',
'dic', 'atmospheric_co2', 'alkalinity', 'hSWS'
])
def set_initial_conditions(self, vs):
# initial conditions for T and S
temp_data = self._read_forcing(vs, 'temperature')[:, :, ::-1]
vs.temp[2:-2, 2:-2, :, :2] = temp_data[:, :, :, np.newaxis] * \
vs.maskT[2:-2, 2:-2, :, np.newaxis]
salt_data = self._read_forcing(vs, 'salinity')[:, :, ::-1]
vs.salt[2:-2, 2:-2, :, :2] = salt_data[..., np.newaxis] * vs.maskT[2:-2, 2:-2, :, np.newaxis]
# use Trenberth wind stress from MITgcm instead of ECMWF (also contained in ecmwf_4deg.cdf)
vs.taux[2:-2, 2:-2, :] = self._read_forcing(vs, 'tau_x')
vs.tauy[2:-2, 2:-2, :] = self._read_forcing(vs, 'tau_y')
# heat flux
with h5netcdf.File(DATA_FILES['ecmwf'], 'r') as ecmwf_data:
qnec_var = ecmwf_data.variables['Q3']
vs.qnec[2:-2, 2:-2, :] = np.array(qnec_var, dtype=str(qnec_var.dtype)).transpose()
vs.qnec[vs.qnec <= -1e10] = 0.0
q = self._read_forcing(vs, 'q_net')
vs.qnet[2:-2, 2:-2, :] = -q
vs.qnet[vs.qnet <= -1e10] = 0.0
fxa = np.sum(vs.qnet[2:-2, 2:-2, :] * vs.area_t[2:-2, 2:-2, np.newaxis]) \
/ 12 / np.sum(vs.area_t[2:-2, 2:-2])
print(' removing an annual mean heat flux imbalance of %e W/m^2' % fxa)
vs.qnet[...] = (vs.qnet - fxa) * vs.maskT[:, :, -1, np.newaxis]
# SST and SSS
vs.sst_clim[2:-2, 2:-2, :] = self._read_forcing(vs, 'sst')
vs.sss_clim[2:-2, 2:-2, :] = self._read_forcing(vs, 'sss')
if vs.enable_idemix:
vs.forc_iw_bottom[2:-2, 2:-2] = self._read_forcing(vs, 'tidal_energy') / vs.rho_0
vs.forc_iw_surface[2:-2, 2:-2] = self._read_forcing(vs, 'wind_energy') / vs.rho_0 * 0.2
if vs.enable_npzd:
phytoplankton = 0.14 * np.exp(vs.zw / 100) * vs.maskT
zooplankton = 0.014 * np.exp(vs.zw / 100) * vs.maskT
vs.phytoplankton[:, :, :, :] = phytoplankton[..., np.newaxis]
vs.zooplankton[:, :, :, :] = zooplankton[..., np.newaxis]
vs.detritus[:, :, :, :] = 1e-4 * vs.maskT[..., np.newaxis]
vs.po4[:, :, :, :] = 2.2
vs.po4[...] *= vs.maskT[..., np.newaxis]
if vs.enable_carbon:
vs.dic[...] = 2300 * vs.maskT[..., np.newaxis]
vs.atmospheric_co2[...] = 280
vs.alkalinity[...] = 2400 * vs.maskT[..., np.newaxis]
vs.hSWS[...] = 5e-7
@veros_method
def set_forcing(self, vs):
year_in_seconds = 360 * 86400.
(n1, f1), (n2, f2) = veros.tools.get_periodic_interval(
vs.time, year_in_seconds, year_in_seconds / 12., 12
)
# wind stress
vs.surface_taux[...] = (f1 * vs.taux[:, :, n1] + f2 * vs.taux[:, :, n2])
vs.surface_tauy[...] = (f1 * vs.tauy[:, :, n1] + f2 * vs.tauy[:, :, n2])
# tke flux
if vs.enable_tke:
vs.forc_tke_surface[1:-1, 1:-1] = np.sqrt((0.5 * (vs.surface_taux[1:-1, 1:-1] \
+ vs.surface_taux[:-2, 1:-1]) / vs.rho_0)**2
+ (0.5 * (vs.surface_tauy[1:-1, 1:-1] \
+ vs.surface_tauy[1:-1, :-2]) / vs.rho_0)**2)**(3. / 2.)
# heat flux : W/m^2 K kg/J m^3/kg = K m/s
cp_0 = 3991.86795711963
sst = f1 * vs.sst_clim[:, :, n1] + f2 * vs.sst_clim[:, :, n2]
qnec = f1 * vs.qnec[:, :, n1] + f2 * vs.qnec[:, :, n2]
qnet = f1 * vs.qnet[:, :, n1] + f2 * vs.qnet[:, :, n2]
vs.forc_temp_surface[...] = (qnet + qnec * (sst - vs.temp[:, :, -1, vs.tau])) \
* vs.maskT[:, :, -1] / cp_0 / vs.rho_0
# salinity restoring
t_rest = 30 * 86400.0
sss = f1 * vs.sss_clim[:, :, n1] + f2 * vs.sss_clim[:, :, n2]
vs.forc_salt_surface[:] = 1. / t_rest * \
(sss - vs.salt[:, :, -1, vs.tau]) * vs.maskT[:, :, -1] * vs.dzt[-1]
# apply simple ice mask
mask = np.logical_and(vs.temp[:, :, -1, vs.tau] * vs.maskT[:, :, -1] < -1.8,
vs.forc_temp_surface < 0.)
vs.forc_temp_surface[mask] = 0.0
vs.forc_salt_surface[mask] = 0.0
if vs.enable_npzd:
# incoming shortwave radiation for plankton production
vs.swr[2:-2, 2:-2] = (f1 * vs.swr_initial[:, :, n1] + f2 * vs.swr_initial[:, :, n2])
@veros_method
def set_diagnostics(self, vs):
vs.diagnostics['snapshot'].output_frequency = 360 * 86400.
vs.diagnostics['overturning'].output_frequency = 360 * 86400.
vs.diagnostics['overturning'].sampling_frequency = vs.dt_tracer
vs.diagnostics['energy'].output_frequency = 360 * 86400.
vs.diagnostics['energy'].sampling_frequency = 86400
average_vars = [
'phytoplankton', 'po4', 'zooplankton', 'detritus', 'wind_speed', 'dic', 'alkalinity',
'temp', 'salt', 'u', 'v', 'surface_tauy', 'surface_taux', 'kappaH', 'psi', 'rho',
'pCO2', 'dpCO2', 'dco2star', 'cflux', 'atmospheric_co2'
]
vs.diagnostics['averages'].output_variables = average_vars
vs.diagnostics['averages'].output_frequency = 360 * 86400.
vs.diagnostics['averages'].sampling_frequency = 86400
vs.diagnostics['npzd'].output_frequency = 10 * 86400
vs.diagnostics['npzd'].save_graph = False
@veros_method
def after_timestep(self, vs):
pass
@veros.tools.cli
def run(*args, **kwargs):
simulation = GlobalFourDegreeBGC(*args, **kwargs)
simulation.setup()
simulation.run()
if __name__ == '__main__':
run()
npzd:
selected_rules:
- "group_npzd_basic"
- "group_carbon_implicit_caco3"
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