Unverified Commit 8496b08a authored by Dion Häfner's avatar Dion Häfner Committed by GitHub

Merge pull request #67 from dionhaefner/potential_density

Potential density is now computed consistently for all EOS
parents c881afd7 d95c01a5
......@@ -23,25 +23,26 @@ def get_rho(vs, salt_loc, temp_loc, press):
@veros_method
def get_potential_rho(vs, salt_loc, temp_loc, press):
def get_potential_rho(vs, salt_loc, temp_loc, press_ref=0.):
"""
calculate potential density as a function of temperature, salinity
and pressure
and reference pressure
Note:
This is identical to get_rho for eq_of_state_type {1, 2, 4}
NB! The potential density is computed for equation of state No.5
according to TEOS-10 formulation, for other equations of state (No.1-4)
it is in-situ density from get_rho
"""
if vs.eq_of_state_type == 1:
return lq.linear_eq_of_state_rho(salt_loc, temp_loc)
elif vs.eq_of_state_type == 2:
return nq1.nonlin1_eq_of_state_rho(salt_loc, temp_loc)
elif vs.eq_of_state_type == 3:
return nq2.nonlin2_eq_of_state_rho(salt_loc, temp_loc, press)
return nq2.nonlin2_eq_of_state_rho(salt_loc, temp_loc, press_ref)
elif vs.eq_of_state_type == 4:
return nq3.nonlin3_eq_of_state_rho(salt_loc, temp_loc)
elif vs.eq_of_state_type == 5:
return gsw.gsw_pot_rho_ct(vs, salt_loc, temp_loc)
return gsw.gsw_rho(vs, salt_loc, temp_loc, press_ref)
else:
raise ValueError('unknown equation of state')
......
This diff is collapsed.
......@@ -174,8 +174,6 @@ def calc_initial_conditions(vs):
vs.rho[...] = density.get_rho(vs, vs.salt, vs.temp, np.abs(vs.zt)[:, np.newaxis]) \
* vs.maskT[..., np.newaxis]
vs.prho[...] = density.get_potential_rho(vs, vs.salt[..., vs.tau], vs.temp[..., vs.tau], np.abs(vs.zt)) \
* vs.maskT[...]
vs.Hd[...] = density.get_dyn_enthalpy(vs, vs.salt, vs.temp, np.abs(vs.zt)[:, np.newaxis]) \
* vs.maskT[..., np.newaxis]
vs.int_drhodT[...] = density.get_int_drhodT(vs, vs.salt, vs.temp, np.abs(vs.zt)[:, np.newaxis])
......
......@@ -253,32 +253,33 @@ def calc_eq_of_state(vs, n):
calculate density, stability frequency, dynamic enthalpy and derivatives
for time level n from temperature and salinity
"""
density_args = (vs, vs.salt[..., n], vs.temp[..., n], np.abs(vs.zt))
salt = vs.salt[..., n]
temp = vs.temp[..., n]
press = np.abs(vs.zt)
"""
calculate new density
"""
vs.rho[..., n] = density.get_rho(*density_args) * vs.maskT
vs.rho[..., n] = density.get_rho(vs, salt, temp, press) * vs.maskT
"""
calculate new potential density
"""
vs.prho[...] = density.get_potential_rho(*density_args) * vs.maskT
vs.prho[...] = density.get_potential_rho(vs, salt, temp) * vs.maskT
"""
calculate new dynamic enthalpy and derivatives
"""
if vs.enable_conserve_energy:
"""
calculate new dynamic enthalpy and derivatives
"""
vs.Hd[..., n] = density.get_dyn_enthalpy(*density_args) * vs.maskT
vs.int_drhodT[..., n] = density.get_int_drhodT(*density_args)
vs.int_drhodS[..., n] = density.get_int_drhodS(*density_args)
vs.Hd[..., n] = density.get_dyn_enthalpy(vs, salt, temp, press) * vs.maskT
vs.int_drhodT[..., n] = density.get_int_drhodT(vs, salt, temp, press)
vs.int_drhodS[..., n] = density.get_int_drhodS(vs, salt, temp, press)
"""
new stability frequency
"""
fxa = -vs.grav / vs.rho_0 / vs.dzw[np.newaxis, np.newaxis, :-1] * vs.maskW[:, :, :-1]
vs.Nsqr[:, :, :-1, n] = fxa * (density.get_rho(
vs, vs.salt[:, :, 1:, n], vs.temp[:, :, 1:, n], np.abs(vs.zt[:-1])
) - vs.rho[:, :, :-1, n]
)
vs, salt[:, :, 1:], temp[:, :, 1:], press[:-1]
) - vs.rho[:, :, :-1, n])
vs.Nsqr[:, :, -1, n] = vs.Nsqr[:, :, -2, n]
......@@ -72,9 +72,9 @@ class Overturning(VerosDiagnostic):
self._allocate(vs)
# sigma levels
self.sige = float(density.get_rho(vs, 35., -2., self.p_ref))
self.sigs = float(density.get_rho(vs, 35., 30., self.p_ref))
self.dsig = (self.sige - self.sigs) / (self.nlevel - 1)
self.sige = density.get_potential_rho(vs, 35., -2., press_ref=self.p_ref)
self.sigs = density.get_potential_rho(vs, 35., 30., press_ref=self.p_ref)
self.dsig = float(self.sige - self.sigs) / (self.nlevel - 1)
logger.debug(' Sigma ranges for overturning diagnostic:')
logger.debug(' Start sigma0 = {:.1f}'.format(self.sigs))
......
......@@ -251,7 +251,7 @@ MAIN_VARIABLES = OrderedDict([
('prho', Variable(
'Potential density', T_GRID, 'kg/m^3',
'Potential density anomaly, relative to the surface mean value of 1024 kg/m^3 '
'(equal to in-situ density anomaly for equation of state 1 to 4)',
'(identical to in-situ density anomaly for equation of state type 1, 2, and 4)',
output=True
)),
......
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