diff --git a/examples/examples.ipynb b/examples/examples.ipynb
index 9b28f1c479444a4c908ae02089c65d957ea3c800..67dd7d382d6fb62580adc0ce4fdab8b24c055d42 100644
--- a/examples/examples.ipynb
+++ b/examples/examples.ipynb
@@ -52,7 +52,7 @@
}
],
"source": [
- "mt.es_liq_hardy(np.asarray([273.16,260.]))"
+ "mt.es_liq_hardy(np.asarray([273.16, 260.0]))"
]
},
{
@@ -82,38 +82,38 @@
}
],
"source": [
- "es = mt.es_liq_analytic\n",
- "p2q = mt.partial_pressure_to_specific_humidity\n",
+ "es = mt.es_liq_analytic\n",
+ "p2q = mt.partial_pressure_to_specific_humidity\n",
"theta_l = mt.theta_l\n",
- "Tl2T = np.vectorize(mt.T_from_Tl)\n",
+ "Tl2T = np.vectorize(mt.T_from_Tl)\n",
"\n",
- "Psfc = 102000.\n",
- "Tsfc = 300.\n",
- "Tmin = 190.\n",
- "dP = 1000.\n",
- "P = np.arange(Psfc,0.4e4,-dP)\n",
+ "Psfc = 102000.0\n",
+ "Tsfc = 300.0\n",
+ "Tmin = 190.0\n",
+ "dP = 1000.0\n",
+ "P = np.arange(Psfc, 0.4e4, -dP)\n",
"\n",
- "RH = 0.77\n",
- "qt = p2q(RH*es(Tsfc),Psfc)\n",
- "print (qt)\n",
+ "RH = 0.77\n",
+ "qt = p2q(RH * es(Tsfc), Psfc)\n",
+ "print(qt)\n",
"\n",
- "sns.set_context('talk')\n",
- "fig, ax = plt.subplots(figsize = (7,5), constrained_layout = True)\n",
+ "sns.set_context(\"talk\")\n",
+ "fig, ax = plt.subplots(figsize=(7, 5), constrained_layout=True)\n",
"\n",
- "Tl = theta_l(Tsfc,Psfc,qt)\n",
- "TK = np.maximum(Tl2T(Tl,P,qt),Tmin)\n",
- "ax.plot(TK,P/100.,label=f\"$\\\\theta_l$ = {Tl:.1f} K, $q_t = ${1000*qt:.2f} g/kg\")\n",
+ "Tl = theta_l(Tsfc, Psfc, qt)\n",
+ "TK = np.maximum(Tl2T(Tl, P, qt), Tmin)\n",
+ "ax.plot(TK, P / 100.0, label=f\"$\\\\theta_l$ = {Tl:.1f} K, $q_t = ${1000*qt:.2f} g/kg\")\n",
"\n",
- "Plcl = mt.plcl(Tsfc,Psfc,qt).squeeze()/100.\n",
+ "Plcl = mt.plcl(Tsfc, Psfc, qt).squeeze() / 100.0\n",
"\n",
- "ax.hlines(Plcl,260,305.,ls='dotted',color='k')\n",
- "ax.set_yticks([Psfc/100,Plcl,850.,700,500,200])\n",
+ "ax.hlines(Plcl, 260, 305.0, ls=\"dotted\", color=\"k\")\n",
+ "ax.set_yticks([Psfc / 100, Plcl, 850.0, 700, 500, 200])\n",
"plt.gca().invert_yaxis()\n",
"plt.legend()\n",
"\n",
"ax.set_xlabel(\"$T$ / K\")\n",
"ax.set_ylabel(\"$P$ / hPa\")\n",
- "#plt.grid()\n",
+ "# plt.grid()\n",
"sns.despine(offset=10)"
]
},
@@ -148,34 +148,38 @@
],
"source": [
"theta_e, theta_l, theta_s = mt.theta_e, mt.theta_l, mt.theta_s\n",
- "Te2T, Tl2T, Ts2T = np.vectorize(mt.T_from_Te), np.vectorize(mt.T_from_Tl), np.vectorize(mt.T_from_Ts)\n",
- "\n",
- "Tmin = 190.\n",
- "dP = 1000.\n",
- "P = np.arange(Psfc,0.4e4,-dP)\n",
- "\n",
- "Psfc = 102000.\n",
- "Tsfc = 300.\n",
- "qt = 15.e-3\n",
- "\n",
- "sns.set_context('talk')\n",
- "fig, ax = plt.subplots(1,2, figsize = (9,7), constrained_layout = True, sharey=True)\n",
- "\n",
- "es = mt.es_liq_analytic\n",
- "TKl = np.maximum(Tl2T(theta_l(Tsfc,Psfc,qt,es=es),P,qt,es=es),Tmin)\n",
- "TKe = np.maximum(Te2T(theta_e(Tsfc,Psfc,qt,es=es),P,qt,es=es),Tmin)\n",
- "TKs = np.maximum(Ts2T(theta_s(Tsfc,Psfc,qt,es=es),P,qt,es=es),Tmin)\n",
- "ax[1].plot(TKe-TKl,P/100.,label=f\"$\\\\theta_e-\\\\theta_l$\")\n",
- "ax[1].plot(TKs-TKl,P/100.,label=f\"$\\\\theta_s-\\\\theta_l$\")\n",
- "ax[1].set_title('es_liq_analytic')\n",
- "\n",
- "es = mt.es_liq\n",
- "TKl = np.maximum(Tl2T(theta_l(Tsfc,Psfc,qt,es=es),P,qt,es=es),Tmin)\n",
- "TKe = np.maximum(Te2T(theta_e(Tsfc,Psfc,qt,es=es),P,qt,es=es),Tmin)\n",
- "TKs = np.maximum(Ts2T(theta_s(Tsfc,Psfc,qt,es=es),P,qt,es=es),Tmin)\n",
- "ax[0].plot(TKe-TKl,P/100.,label=f\"$\\\\theta_e-\\\\theta_l$\")\n",
- "ax[0].plot(TKs-TKl,P/100.,label=f\"$\\\\theta_s-\\\\theta_l$\")\n",
- "ax[0].set_title('es_liq')\n",
+ "Te2T, Tl2T, Ts2T = (\n",
+ " np.vectorize(mt.T_from_Te),\n",
+ " np.vectorize(mt.T_from_Tl),\n",
+ " np.vectorize(mt.T_from_Ts),\n",
+ ")\n",
+ "\n",
+ "Tmin = 190.0\n",
+ "dP = 1000.0\n",
+ "P = np.arange(Psfc, 0.4e4, -dP)\n",
+ "\n",
+ "Psfc = 102000.0\n",
+ "Tsfc = 300.0\n",
+ "qt = 15.0e-3\n",
+ "\n",
+ "sns.set_context(\"talk\")\n",
+ "fig, ax = plt.subplots(1, 2, figsize=(9, 7), constrained_layout=True, sharey=True)\n",
+ "\n",
+ "es = mt.es_liq_analytic\n",
+ "TKl = np.maximum(Tl2T(theta_l(Tsfc, Psfc, qt, es=es), P, qt, es=es), Tmin)\n",
+ "TKe = np.maximum(Te2T(theta_e(Tsfc, Psfc, qt, es=es), P, qt, es=es), Tmin)\n",
+ "TKs = np.maximum(Ts2T(theta_s(Tsfc, Psfc, qt, es=es), P, qt, es=es), Tmin)\n",
+ "ax[1].plot(TKe - TKl, P / 100.0, label=f\"$\\\\theta_e-\\\\theta_l$\")\n",
+ "ax[1].plot(TKs - TKl, P / 100.0, label=f\"$\\\\theta_s-\\\\theta_l$\")\n",
+ "ax[1].set_title(\"es_liq_analytic\")\n",
+ "\n",
+ "es = mt.es_liq\n",
+ "TKl = np.maximum(Tl2T(theta_l(Tsfc, Psfc, qt, es=es), P, qt, es=es), Tmin)\n",
+ "TKe = np.maximum(Te2T(theta_e(Tsfc, Psfc, qt, es=es), P, qt, es=es), Tmin)\n",
+ "TKs = np.maximum(Ts2T(theta_s(Tsfc, Psfc, qt, es=es), P, qt, es=es), Tmin)\n",
+ "ax[0].plot(TKe - TKl, P / 100.0, label=f\"$\\\\theta_e-\\\\theta_l$\")\n",
+ "ax[0].plot(TKs - TKl, P / 100.0, label=f\"$\\\\theta_s-\\\\theta_l$\")\n",
+ "ax[0].set_title(\"es_liq\")\n",
"\n",
"plt.gca().invert_yaxis()\n",
"\n",
@@ -206,9 +210,9 @@
"outputs": [],
"source": [
"# Version 1.0 released by David Romps on September 12, 2017.\n",
- "# \n",
+ "#\n",
"# When using this code, please cite:\n",
- "# \n",
+ "#\n",
"# @article{16lcl,\n",
"# Title = {Exact expression for the lifting condensation level},\n",
"# Author = {David M. Romps},\n",
@@ -224,46 +228,56 @@
"# - Exactly one of rh, rhl, and rhs (dimensionless, from 0 to 1):\n",
"# * The value of rh is interpreted to be the relative humidity with\n",
"# respect to liquid water if T >= 273.15 K and with respect to ice if\n",
- "# T < 273.15 K. \n",
+ "# T < 273.15 K.\n",
"# * The value of rhl is interpreted to be the relative humidity with\n",
"# respect to liquid water\n",
"# * The value of rhs is interpreted to be the relative humidity with\n",
"# respect to ice\n",
"# - ldl is an optional logical flag. If true, the lifting deposition\n",
- "# level (LDL) is returned instead of the LCL. \n",
+ "# level (LDL) is returned instead of the LCL.\n",
"# - min_lcl_ldl is an optional logical flag. If true, the minimum of the\n",
"# LCL and LDL is returned.\n",
"\n",
- "def lcl(p,T,rh=None,rhl=None,rhs=None,return_ldl=False,return_min_lcl_ldl=False):\n",
+ "\n",
+ "def lcl(p, T, rh=None, rhl=None, rhs=None, return_ldl=False, return_min_lcl_ldl=False):\n",
"\n",
" import math\n",
" import scipy.special\n",
"\n",
" # Parameters\n",
- " Ttrip = 273.16 # K\n",
- " ptrip = 611.65 # Pa\n",
- " E0v = 2.3740e6 # J/kg\n",
- " E0s = 0.3337e6 # J/kg\n",
- " ggr = 9.81 # m/s^2\n",
- " rgasa = 287.04 # J/kg/K \n",
- " rgasv = 461 # J/kg/K \n",
- " cva = 719 # J/kg/K\n",
- " cvv = 1418 # J/kg/K \n",
- " cvl = 4119 # J/kg/K \n",
- " cvs = 1861 # J/kg/K \n",
- " cpa = cva + rgasa\n",
- " cpv = cvv + rgasv\n",
+ " Ttrip = 273.16 # K\n",
+ " ptrip = 611.65 # Pa\n",
+ " E0v = 2.3740e6 # J/kg\n",
+ " E0s = 0.3337e6 # J/kg\n",
+ " ggr = 9.81 # m/s^2\n",
+ " rgasa = 287.04 # J/kg/K\n",
+ " rgasv = 461 # J/kg/K\n",
+ " cva = 719 # J/kg/K\n",
+ " cvv = 1418 # J/kg/K\n",
+ " cvl = 4119 # J/kg/K\n",
+ " cvs = 1861 # J/kg/K\n",
+ " cpa = cva + rgasa\n",
+ " cpv = cvv + rgasv\n",
"\n",
" # The saturation vapor pressure over liquid water\n",
" def pvstarl(T):\n",
- " return ptrip * (T/Ttrip)**((cpv-cvl)/rgasv) * math.exp( (E0v - (cvv-cvl)*Ttrip) / rgasv * (1/Ttrip - 1/T) )\n",
+ " return (\n",
+ " ptrip\n",
+ " * (T / Ttrip) ** ((cpv - cvl) / rgasv)\n",
+ " * math.exp((E0v - (cvv - cvl) * Ttrip) / rgasv * (1 / Ttrip - 1 / T))\n",
+ " )\n",
+ "\n",
" # The saturation vapor pressure over solid ice\n",
" def pvstars(T):\n",
- " return ptrip * (T/Ttrip)**((cpv-cvs)/rgasv) * math.exp( (E0v + E0s - (cvv-cvs)*Ttrip) / rgasv * (1/Ttrip - 1/T)) \n",
+ " return (\n",
+ " ptrip\n",
+ " * (T / Ttrip) ** ((cpv - cvs) / rgasv)\n",
+ " * math.exp((E0v + E0s - (cvv - cvs) * Ttrip) / rgasv * (1 / Ttrip - 1 / T))\n",
+ " )\n",
"\n",
" # Calculate pv from rh, rhl, or rhs\n",
" rh_counter = 0\n",
- " if rh is not None:\n",
+ " if rh is not None:\n",
" rh_counter = rh_counter + 1\n",
" if rhl is not None:\n",
" rh_counter = rh_counter + 1\n",
@@ -271,11 +285,11 @@
" rh_counter = rh_counter + 1\n",
" if rh_counter != 1:\n",
" print(rh_counter)\n",
- " exit('Error in lcl: Exactly one of rh, rhl, and rhs must be specified')\n",
+ " exit(\"Error in lcl: Exactly one of rh, rhl, and rhs must be specified\")\n",
" if rh is not None:\n",
- " # The variable rh is assumed to be \n",
- " # with respect to liquid if T > Ttrip and \n",
- " # with respect to solid if T < Ttrip\n",
+ " # The variable rh is assumed to be\n",
+ " # with respect to liquid if T > Ttrip and\n",
+ " # with respect to solid if T < Ttrip\n",
" if T > Ttrip:\n",
" pv = rh * pvstarl(T)\n",
" else:\n",
@@ -299,34 +313,34 @@
" if pv > p:\n",
" return N\n",
"\n",
- "# Calculate lcl_liquid and lcl_solid\n",
- " qv = rgasa*pv / (rgasv*p + (rgasa-rgasv)*pv)\n",
- " rgasm = (1-qv)*rgasa + qv*rgasv\n",
- " cpm = (1-qv)*cpa + qv*cpv\n",
+ " # Calculate lcl_liquid and lcl_solid\n",
+ " qv = rgasa * pv / (rgasv * p + (rgasa - rgasv) * pv)\n",
+ " rgasm = (1 - qv) * rgasa + qv * rgasv\n",
+ " cpm = (1 - qv) * cpa + qv * cpv\n",
" if rh == 0:\n",
- " return cpm*T/ggr\n",
- " aL = -(cpv-cvl)/rgasv + cpm/rgasm\n",
- " bL = -(E0v-(cvv-cvl)*Ttrip)/(rgasv*T)\n",
- " cL = pv/pvstarl(T)*math.exp(-(E0v-(cvv-cvl)*Ttrip)/(rgasv*T))\n",
- " aS = -(cpv-cvs)/rgasv + cpm/rgasm\n",
- " bS = -(E0v+E0s-(cvv-cvs)*Ttrip)/(rgasv*T)\n",
- " cS = pv/pvstars(T)*math.exp(-(E0v+E0s-(cvv-cvs)*Ttrip)/(rgasv*T))\n",
- " X = bL/(aL*scipy.special.lambertw(bL/aL*cL**(1/aL),-1).real)\n",
- " Y = bS/(aS*scipy.special.lambertw(bS/aS*cS**(1/aS),-1).real) \n",
- " \n",
- " lcl = cpm*T/ggr*( 1 - X)\n",
- " ldl = cpm*T/ggr*( 1 - Y)\n",
+ " return cpm * T / ggr\n",
+ " aL = -(cpv - cvl) / rgasv + cpm / rgasm\n",
+ " bL = -(E0v - (cvv - cvl) * Ttrip) / (rgasv * T)\n",
+ " cL = pv / pvstarl(T) * math.exp(-(E0v - (cvv - cvl) * Ttrip) / (rgasv * T))\n",
+ " aS = -(cpv - cvs) / rgasv + cpm / rgasm\n",
+ " bS = -(E0v + E0s - (cvv - cvs) * Ttrip) / (rgasv * T)\n",
+ " cS = pv / pvstars(T) * math.exp(-(E0v + E0s - (cvv - cvs) * Ttrip) / (rgasv * T))\n",
+ " X = bL / (aL * scipy.special.lambertw(bL / aL * cL ** (1 / aL), -1).real)\n",
+ " Y = bS / (aS * scipy.special.lambertw(bS / aS * cS ** (1 / aS), -1).real)\n",
+ "\n",
+ " lcl = cpm * T / ggr * (1 - X)\n",
+ " ldl = cpm * T / ggr * (1 - Y)\n",
"\n",
" # Modifications of the code to output Plcl or Pldl\n",
- " Plcl = PPa * X**(cpm/rgasm)\n",
- " Pldl = PPa * X**(cpm/rgasm)\n",
+ " Plcl = PPa * X ** (cpm / rgasm)\n",
+ " Pldl = PPa * X ** (cpm / rgasm)\n",
" # Return either lcl or ldl\n",
" if return_ldl and return_min_lcl_ldl:\n",
- " exit('return_ldl and return_min_lcl_ldl cannot both be true')\n",
+ " exit(\"return_ldl and return_min_lcl_ldl cannot both be true\")\n",
" elif return_ldl:\n",
" return Pldl\n",
" elif return_min_lcl_ldl:\n",
- " return min(Plcl,Pldl)\n",
+ " return min(Plcl, Pldl)\n",
" else:\n",
" return Plcl"
]
@@ -360,54 +374,54 @@
}
],
"source": [
- "PPa = 101325.\n",
- "qt = np.arange(2.5e-3,8.e-3,0.2e-3)\n",
- "TK = 285.\n",
+ "PPa = 101325.0\n",
+ "qt = np.arange(2.5e-3, 8.0e-3, 0.2e-3)\n",
+ "TK = 285.0\n",
"Plcl_X = np.zeros(len(qt))\n",
"Plcl_B = np.zeros(len(qt))\n",
"Plcl_R = np.zeros(len(qt))\n",
"\n",
- "for i,x in enumerate(qt):\n",
- " RH = mt.mixing_ratio_to_partial_pressure(x/(1.-x),PPa)/mt.es_liq(TK)\n",
- " Plcl_R[i] = lcl(PPa,TK,RH)\n",
- " Plcl_X[i] = mt.plcl(TK,PPa,x)\n",
- " Plcl_B[i] = mt.plcl_bolton(TK,PPa,x)\n",
+ "for i, x in enumerate(qt):\n",
+ " RH = mt.mixing_ratio_to_partial_pressure(x / (1.0 - x), PPa) / mt.es_liq(TK)\n",
+ " Plcl_R[i] = lcl(PPa, TK, RH)\n",
+ " Plcl_X[i] = mt.plcl(TK, PPa, x)\n",
+ " Plcl_B[i] = mt.plcl_bolton(TK, PPa, x)\n",
"\n",
- "del1 = (Plcl_B-Plcl_X)/100.\n",
- "del2 = (Plcl_R-Plcl_X)/100.\n",
+ "del1 = (Plcl_B - Plcl_X) / 100.0\n",
+ "del2 = (Plcl_R - Plcl_X) / 100.0\n",
"\n",
- "fig, axs = plt.subplots(1,2, figsize = (7,5), constrained_layout = True, sharey=True)\n",
+ "fig, axs = plt.subplots(1, 2, figsize=(7, 5), constrained_layout=True, sharey=True)\n",
"\n",
- "axs[0].plot(qt*1.e3,del1,label='$\\\\delta_\\mathrm{B}$')\n",
- "axs[0].plot(qt*1.e3,del2,label='$\\\\delta_\\mathrm{R}$')\n",
+ "axs[0].plot(qt * 1.0e3, del1, label=\"$\\\\delta_\\mathrm{B}$\")\n",
+ "axs[0].plot(qt * 1.0e3, del2, label=\"$\\\\delta_\\mathrm{R}$\")\n",
"axs[0].legend(loc=\"best\")\n",
- "axs[0].set_ylabel('$\\delta P$ / hPa')\n",
- "axs[0].set_xlabel('$q_\\mathrm{t}$ / gkg$^{-1}$')\n",
- "axs[0].set_title(f'T={TK:.0f} K')\n",
+ "axs[0].set_ylabel(\"$\\delta P$ / hPa\")\n",
+ "axs[0].set_xlabel(\"$q_\\mathrm{t}$ / gkg$^{-1}$\")\n",
+ "axs[0].set_title(f\"T={TK:.0f} K\")\n",
"\n",
- "qt = np.arange(0.5e-3,28.e-3,0.2e-3)\n",
- "TK = 310.\n",
+ "qt = np.arange(0.5e-3, 28.0e-3, 0.2e-3)\n",
+ "TK = 310.0\n",
"\n",
"Plcl_X = np.zeros(len(qt))\n",
"Plcl_B = np.zeros(len(qt))\n",
"Plcl_R = np.zeros(len(qt))\n",
- "for i,x in enumerate(qt):\n",
- " RH = mt.mixing_ratio_to_partial_pressure(x/(1.-x),PPa)/mt.es_liq(TK)\n",
- " Plcl_R[i] = lcl(PPa,TK,RH)\n",
- " Plcl_X[i] = mt.plcl(TK,PPa,x)\n",
- " Plcl_B[i] = mt.plcl_bolton(TK,PPa,x)\n",
+ "for i, x in enumerate(qt):\n",
+ " RH = mt.mixing_ratio_to_partial_pressure(x / (1.0 - x), PPa) / mt.es_liq(TK)\n",
+ " Plcl_R[i] = lcl(PPa, TK, RH)\n",
+ " Plcl_X[i] = mt.plcl(TK, PPa, x)\n",
+ " Plcl_B[i] = mt.plcl_bolton(TK, PPa, x)\n",
"\n",
- "del1 = (Plcl_B-Plcl_X)/100.\n",
- "del2 = (Plcl_R-Plcl_X)/100.\n",
+ "del1 = (Plcl_B - Plcl_X) / 100.0\n",
+ "del2 = (Plcl_R - Plcl_X) / 100.0\n",
"\n",
- "axs[1].plot(qt*1.e3,del1)\n",
- "axs[1].plot(qt*1.e3,del2)\n",
- "axs[1].set_xlabel('$q_\\mathrm{t}$ / gkg$^{-1}$')\n",
- "axs[1].set_title(f'T={TK:.0f} K')\n",
+ "axs[1].plot(qt * 1.0e3, del1)\n",
+ "axs[1].plot(qt * 1.0e3, del2)\n",
+ "axs[1].set_xlabel(\"$q_\\mathrm{t}$ / gkg$^{-1}$\")\n",
+ "axs[1].set_title(f\"T={TK:.0f} K\")\n",
"\n",
"sns.set_context(\"talk\", font_scale=1.2)\n",
"sns.despine(offset=10)\n",
- "#fig.savefig(plot_dir+'Plcl.pdf')"
+ "# fig.savefig(plot_dir+'Plcl.pdf')"
]
},
{
@@ -451,25 +465,43 @@
}
],
"source": [
- "qt = 16.774409538883497e-3\n",
- "Tice,Py= mt.moist_adiabat(Tsfc,Psfc,150e2,-10.,qt,cc=constants.ci,l=mt.sublimation_enthalpy ,es = mt.es_mxd) \n",
- "Tliq,Px= mt.moist_adiabat(Tsfc,Psfc,150e2,-10.,qt,cc=constants.cl,l=mt.vaporization_enthalpy,es = mt.es_mxd)\n",
- "\n",
- "Tx = np.ones(len(Px))*constants.T0\n",
- "Tx[Tliq>constants.T0] = Tliq[Tliq>constants.T0]\n",
- "Tx[Tice<constants.T0] = Tice[Tice<constants.T0]\n",
- "Tx = np.maximum(Tx,190.)\n",
- "Tliq = np.maximum(Tliq,190.)\n",
- "\n",
- "Tl2T = np.vectorize(mt.T_from_Tl)\n",
- "Tl = mt.theta_l(Tsfc,Psfc,qt)\n",
- "TK = np.maximum(Tl2T(Tl,Px,qt),Tmin)\n",
- "\n",
- "sns.set_context('talk')\n",
- "fig, ax = plt.subplots(figsize = (4,5), constrained_layout = True, sharey=True)\n",
- "\n",
- "ax.plot(mt.theta(Tliq,Px) ,Px/100.,label=f\"liquid\",c='dodgerblue')\n",
- "ax.plot(mt.theta(Tx,Px) ,Px/100.,label=f\"with ice\",c='k',ls='dotted')\n",
+ "qt = 16.774409538883497e-3\n",
+ "Tice, Py = mt.moist_adiabat(\n",
+ " Tsfc,\n",
+ " Psfc,\n",
+ " 150e2,\n",
+ " -10.0,\n",
+ " qt,\n",
+ " cc=constants.ci,\n",
+ " l=mt.sublimation_enthalpy,\n",
+ " es=mt.es_mxd,\n",
+ ")\n",
+ "Tliq, Px = mt.moist_adiabat(\n",
+ " Tsfc,\n",
+ " Psfc,\n",
+ " 150e2,\n",
+ " -10.0,\n",
+ " qt,\n",
+ " cc=constants.cl,\n",
+ " l=mt.vaporization_enthalpy,\n",
+ " es=mt.es_mxd,\n",
+ ")\n",
+ "\n",
+ "Tx = np.ones(len(Px)) * constants.T0\n",
+ "Tx[Tliq > constants.T0] = Tliq[Tliq > constants.T0]\n",
+ "Tx[Tice < constants.T0] = Tice[Tice < constants.T0]\n",
+ "Tx = np.maximum(Tx, 190.0)\n",
+ "Tliq = np.maximum(Tliq, 190.0)\n",
+ "\n",
+ "Tl2T = np.vectorize(mt.T_from_Tl)\n",
+ "Tl = mt.theta_l(Tsfc, Psfc, qt)\n",
+ "TK = np.maximum(Tl2T(Tl, Px, qt), Tmin)\n",
+ "\n",
+ "sns.set_context(\"talk\")\n",
+ "fig, ax = plt.subplots(figsize=(4, 5), constrained_layout=True, sharey=True)\n",
+ "\n",
+ "ax.plot(mt.theta(Tliq, Px), Px / 100.0, label=f\"liquid\", c=\"dodgerblue\")\n",
+ "ax.plot(mt.theta(Tx, Px), Px / 100.0, label=f\"with ice\", c=\"k\", ls=\"dotted\")\n",
"\n",
"plt.gca().invert_yaxis()\n",
"\n",
@@ -495,8 +527,8 @@
}
],
"source": [
- "q=[0.,0.,0]\n",
- "print (np.sum(q))"
+ "q = [0.0, 0.0, 0]\n",
+ "print(np.sum(q))"
]
},
{
diff --git a/moist_thermodynamics/constants.py b/moist_thermodynamics/constants.py
index 023308f56181cb076534c6b8eabb7866fb905498..0055d4a6d91e1d86c8a57bcd0efdc57bc1133758 100644
--- a/moist_thermodynamics/constants.py
+++ b/moist_thermodynamics/constants.py
@@ -5,78 +5,98 @@ Author: Bjorn Stevens (bjorn.stevens@mpimet.mpg.de)
#
import numpy as np
-c = speed_of_light = 299792458
-kB = boltzmann_constant = 1.380649e-23
-N_avo = avogadro_number = 6.02214076e+23
-h = planck_constant = 6.62607015e-34
-G = gravitational_constant = 6.67430e-11
-
-mass_earth = 5.9722e+24
-area_earth = 510065623e6
-gravity_earth = 9.80665 # 4*np.pi*G*mass_earth/area_earth gives a value of 9.8203
-radius_earth = np.sqrt(G*mass_earth/gravity_earth) # 6375.4, where standard equatorial radius is given as 6,378.1 km
-
-Rstar = ideal_gas_constant = kB*N_avo
-sigma = stefan_boltzmann_constant = 2*(np.pi**5) * kB**4 / (15 * h**3 * c**2)
-
-Tsfc = mean_earth_surface_temperature = 288. # Kelvin
-Psfc = mean_earth_surface_pressure = 98443 # Pa
-mass_atm = mass_earth_atmosphere = (Psfc / gravity_earth) # kg/m^2
-
-P0 = standard_pressure = 100000. # Pa
-T0 = standard_temperature = 273.15 # K
-P00 = standard_atmosphere_pressure = 101325. # Pa
-T00 = standard_atmosphere_temperature = 298.15 # K
-Ttriple = temperature_triple_point = 273.16 # K
+c = speed_of_light = 299792458
+kB = boltzmann_constant = 1.380649e-23
+N_avo = avogadro_number = 6.02214076e23
+h = planck_constant = 6.62607015e-34
+G = gravitational_constant = 6.67430e-11
+
+mass_earth = 5.9722e24
+area_earth = 510065623e6
+gravity_earth = 9.80665 # 4*np.pi*G*mass_earth/area_earth gives a value of 9.8203
+radius_earth = np.sqrt(
+ G * mass_earth / gravity_earth
+) # 6375.4, where standard equatorial radius is given as 6,378.1 km
+
+Rstar = ideal_gas_constant = kB * N_avo
+sigma = stefan_boltzmann_constant = 2 * (np.pi**5) * kB**4 / (15 * h**3 * c**2)
+
+Tsfc = mean_earth_surface_temperature = 288.0 # Kelvin
+Psfc = mean_earth_surface_pressure = 98443 # Pa
+mass_atm = mass_earth_atmosphere = Psfc / gravity_earth # kg/m^2
+
+P0 = standard_pressure = 100000.0 # Pa
+T0 = standard_temperature = 273.15 # K
+P00 = standard_atmosphere_pressure = 101325.0 # Pa
+T00 = standard_atmosphere_temperature = 298.15 # K
+Ttriple = temperature_triple_point = 273.16 # K
#
# Based on Park et al (2004) Meteorlogia, O2 levels are declining as CO2 levels rise, but at a tiny rate.
#
-x_ar = earth_atmosphere_ar_mass_fraction = 9.332e-3
-x_o2 = earth_atmosphere_o2_mass_fraction = 0.20944
-x_n2 = earth_atmosphere_n2_mass_fraction = 0.78083
+x_ar = earth_atmosphere_ar_mass_fraction = 9.332e-3
+x_o2 = earth_atmosphere_o2_mass_fraction = 0.20944
+x_n2 = earth_atmosphere_n2_mass_fraction = 0.78083
-x_co2 = earth_atmosphere_co2_mass_fraction = 0.415e-3
-x_ch4 = earth_atmosphere_methane_mass_fraction = 772.e-9
-x_n2o = earth_atmosphere_nitrous_oxide_mass_fraction = 334.e-9
-x_o3 = earth_atmosphere_ozone_mass_fraction = 200.e-9
+x_co2 = earth_atmosphere_co2_mass_fraction = 0.415e-3
+x_ch4 = earth_atmosphere_methane_mass_fraction = 772.0e-9
+x_n2o = earth_atmosphere_nitrous_oxide_mass_fraction = 334.0e-9
+x_o3 = earth_atmosphere_ozone_mass_fraction = 200.0e-9
#
# Based on Chase (1998) J Phys Chem Ref Data
#
-m_ar = molar_mass_ar = 39.948
-m_o2 = molar_mass_o2 = 15.9994 * 2
-m_n2 = molar_mass_n2 = 14.0067 * 2
+m_ar = molar_mass_ar = 39.948
+m_o2 = molar_mass_o2 = 15.9994 * 2
+m_n2 = molar_mass_n2 = 14.0067 * 2
m_co2 = molar_mass_co2 = 44.011
m_h2o = molar_mass_h2o = 18.01528
-m_ch4 = molar_mass_methane = m_co2 + 4*m_h2o - 3.*m_o2
-m_n2o = molar_mass_nitroush_oxide = m_n2 + m_o2/2
-m_o3 = molar_mass_nitroush_ozone = m_o2 * 3/2
+m_ch4 = molar_mass_methane = m_co2 + 4 * m_h2o - 3.0 * m_o2
+m_n2o = molar_mass_nitroush_oxide = m_n2 + m_o2 / 2
+m_o3 = molar_mass_nitroush_ozone = m_o2 * 3 / 2
# molar mass of dry air
-md = atomic_mass_dry_air = x_ar*m_ar + x_o2*m_o2 + x_n2*m_n2 + x_co2*m_co2 + x_ch4*m_ch4 + x_n2o*m_n2o + x_o3*m_o3
-
-cp_ar = isobaric_specific_heat_capacity_a = 20.786 # 298.15K
-cp_o2 = isobaric_specific_heat_capacity_o2 = 29.376 # isobaric_specific_heat_capacity_oxygen = 298.15K or 29.126 @ 200K
-cp_n2 = isobaric_specific_heat_capacity_n2 = 29.124 # 298.15K or 29.107 @ 200K
-cp_co2 = isobaric_specific_heat_capacity_co2= 37.129 # 298.15K or 32.359 @ 200K
-cp_h2o = isobaric_specific_heat_capacity_h2o= 33.349 + (33.590 - 33.349)/98.15 * (T0-200) # Interpolated to T0 from Chase values (but not used)
-
-s0_ar = entropy_argon_stp = 154.845 # J/mol*K
-s0_o2 = entropy_oxygen_stp = 205.147 # J/mol*K
-s0_n2 = entropy_nitrogen_stp = 191.609 # J/mol*K
-s0_co2 = entropy_co2_stp = 213.795 # J/mol*K
-s0_h2o = entropy_h2o_stp = 188.854 # J/mol*K
-
-q_ar = argon_mass_mixing_fraction = x_ar *m_ar /md
-q_o2 = o2_mass_mixing_fraction = x_o2 *m_o2 /md
-q_n2 = n2_mass_mixing_fraction = x_n2 *m_n2 /md
-q_co2 = water_vapor_mixing_fraction = x_co2*m_co2/md
-
-Rd = dry_air_gas_constant = (Rstar/md)*(x_ar+x_o2+x_n2+x_co2) * 1000. # J/kg/K
-cpd = isobaric_dry_air_specific_heat = ( 1./md)*(x_ar*cp_ar + x_o2*cp_o2 + x_n2*cp_n2 + x_co2*cp_co2) *1000. # J/kg/K
-sd0 = entropy_dry_air_stp = ( 1./md)*(x_ar*s0_ar + x_o2*s0_o2 + x_n2*s0_n2 + x_co2*s0_co2) *1000. # J/kg*K
-sd00 = entropy_dry_air_satmt = sd0 + cpd * np.log(T0/T00)
+md = atomic_mass_dry_air = (
+ x_ar * m_ar
+ + x_o2 * m_o2
+ + x_n2 * m_n2
+ + x_co2 * m_co2
+ + x_ch4 * m_ch4
+ + x_n2o * m_n2o
+ + x_o3 * m_o3
+)
+
+cp_ar = isobaric_specific_heat_capacity_a = 20.786 # 298.15K
+cp_o2 = (
+ isobaric_specific_heat_capacity_o2
+) = 29.376 # isobaric_specific_heat_capacity_oxygen = 298.15K or 29.126 @ 200K
+cp_n2 = isobaric_specific_heat_capacity_n2 = 29.124 # 298.15K or 29.107 @ 200K
+cp_co2 = isobaric_specific_heat_capacity_co2 = 37.129 # 298.15K or 32.359 @ 200K
+cp_h2o = isobaric_specific_heat_capacity_h2o = 33.349 + (33.590 - 33.349) / 98.15 * (
+ T0 - 200
+) # Interpolated to T0 from Chase values (but not used)
+
+s0_ar = entropy_argon_stp = 154.845 # J/mol*K
+s0_o2 = entropy_oxygen_stp = 205.147 # J/mol*K
+s0_n2 = entropy_nitrogen_stp = 191.609 # J/mol*K
+s0_co2 = entropy_co2_stp = 213.795 # J/mol*K
+s0_h2o = entropy_h2o_stp = 188.854 # J/mol*K
+
+q_ar = argon_mass_mixing_fraction = x_ar * m_ar / md
+q_o2 = o2_mass_mixing_fraction = x_o2 * m_o2 / md
+q_n2 = n2_mass_mixing_fraction = x_n2 * m_n2 / md
+q_co2 = water_vapor_mixing_fraction = x_co2 * m_co2 / md
+
+Rd = dry_air_gas_constant = (
+ (Rstar / md) * (x_ar + x_o2 + x_n2 + x_co2) * 1000.0
+) # J/kg/K
+cpd = isobaric_dry_air_specific_heat = (
+ (1.0 / md) * (x_ar * cp_ar + x_o2 * cp_o2 + x_n2 * cp_n2 + x_co2 * cp_co2) * 1000.0
+) # J/kg/K
+sd0 = entropy_dry_air_stp = (
+ (1.0 / md) * (x_ar * s0_ar + x_o2 * s0_o2 + x_n2 * s0_n2 + x_co2 * s0_co2) * 1000.0
+) # J/kg*K
+sd00 = entropy_dry_air_satmt = sd0 + cpd * np.log(T0 / T00)
#
# cl and ci, especially ci, varies considerably with temperature. Consider that
# cl = 4273 J/kg/K at 263 K decreases sharply to 4220 J/kg/K by 273 K and ever more slowly to
@@ -86,27 +106,37 @@ sd00 = entropy_dry_air_satmt = sd0 + cpd * np.log(T0/T00)
# At standard temperature and pressure they have the values
# cl = 4219.32 # ''
# ci = 2096.70 # ''
-cpv = isobaric_water_vapor_specific_heat = 1865.01 # IAPWS97 at 273.15
-cl = liquid_water_specific_heat = 4179.57 # IAPWS97 at 305 and P=0.1 MPa (gives a good fit for es over ice)
-ci = frozen_water_specific_heat = 1905.43 # IAPWS97 at 247.065 and P=0.1 MPa (gives a good fit for es over ice)
-delta_cl= cpv-cl
-delta_ci= cpv-ci
-
-lv0 = vaporization_enthalpy_stp = 2500.93e3 # IAPWS97 at 273.15
-lf0 = melting_enthalpy_stp = 333.42e3 # ''
-ls0 = sublimation_enthalpy_stp = lv0+lf0 # ''
-
-Rv = water_vapor_gas_constant = (Rstar/m_h2o) *1000. #J/kg/K
-sv00 = entropy_water_vapor_satmt = (s0_h2o/m_h2o)*1000. + cpv * np.log(T0/298.15)
-
-eps1 = rd_over_rv = Rd/Rv
-eps2 = rv_over_rd_minus_one = Rv/Rd -1.
-
-TvC = temperature_water_vapor_critical_point = 647.096 # Critical temperature [K] of water vapor
-PvC = pressure_water_vapor_critical_point = 22.064e6 # Critical pressure [Pa] of water vapor
-
-TvT = temperature_water_vapor_triple_point = 273.16 # Triple point temperature [K] of water
-PvT = pressure_water_vapor_triple_point = 611.655
-lvT = vaporization_enthalpy_triple_point = lv0 + (cpv-cl)*(TvT-T0)
-lfT = melting_enthalpy_triple_point = lf0 + (cpv-ci)*(TvT-T0)
-lsT = sublimation_enthalpy_triple_point = lvT + lfT
+cpv = isobaric_water_vapor_specific_heat = 1865.01 # IAPWS97 at 273.15
+cl = (
+ liquid_water_specific_heat
+) = 4179.57 # IAPWS97 at 305 and P=0.1 MPa (gives a good fit for es over ice)
+ci = (
+ frozen_water_specific_heat
+) = 1905.43 # IAPWS97 at 247.065 and P=0.1 MPa (gives a good fit for es over ice)
+delta_cl = cpv - cl
+delta_ci = cpv - ci
+
+lv0 = vaporization_enthalpy_stp = 2500.93e3 # IAPWS97 at 273.15
+lf0 = melting_enthalpy_stp = 333.42e3 # ''
+ls0 = sublimation_enthalpy_stp = lv0 + lf0 # ''
+
+Rv = water_vapor_gas_constant = (Rstar / m_h2o) * 1000.0 # J/kg/K
+sv00 = entropy_water_vapor_satmt = (s0_h2o / m_h2o) * 1000.0 + cpv * np.log(T0 / 298.15)
+
+eps1 = rd_over_rv = Rd / Rv
+eps2 = rv_over_rd_minus_one = Rv / Rd - 1.0
+
+TvC = (
+ temperature_water_vapor_critical_point
+) = 647.096 # Critical temperature [K] of water vapor
+PvC = (
+ pressure_water_vapor_critical_point
+) = 22.064e6 # Critical pressure [Pa] of water vapor
+
+TvT = (
+ temperature_water_vapor_triple_point
+) = 273.16 # Triple point temperature [K] of water
+PvT = pressure_water_vapor_triple_point = 611.655
+lvT = vaporization_enthalpy_triple_point = lv0 + (cpv - cl) * (TvT - T0)
+lfT = melting_enthalpy_triple_point = lf0 + (cpv - ci) * (TvT - T0)
+lsT = sublimation_enthalpy_triple_point = lvT + lfT
diff --git a/moist_thermodynamics/functions.py b/moist_thermodynamics/functions.py
index 3ba9178a6878dbc0c521ca449659e854a4fd6e96..42cf33067f84e87f6196f9ec47ec766baed91968 100644
--- a/moist_thermodynamics/functions.py
+++ b/moist_thermodynamics/functions.py
@@ -12,60 +12,73 @@ from . import constants
import numpy as np
from scipy import interpolate, optimize
-def planck(T,nu):
+
+def planck(T, nu):
"""Planck source function (J/m2 per steradian per Hz)
-
+
Args:
T: temperature in kelvin
nu: frequency in Hz
-
+
Returns:
Returns the radiance in the differential frequency interval per unit steradian. Usually we
- multiply by $\pi$ to convert to irradiances
-
+ multiply by $\pi$ to convert to irradiances
+
>>> planck(300,1000*constants.c)
8.086837160291128e-15
"""
c = constants.speed_of_light
h = constants.planck_constant
- kB= constants.boltzmann_constant
- return (2 * h * nu**3 / c**2) / (np.exp(h*nu/(kB*T))-1)
+ kB = constants.boltzmann_constant
+ return (2 * h * nu**3 / c**2) / (np.exp(h * nu / (kB * T)) - 1)
+
def es_liq(T):
"""Returns saturation vapor pressure (Pa) over planer liquid water
-
+
Encodes the empirical fits of Wagner and Pruss (2002), Eq 2.5a (page 399). Their formulation
- is compared to other fits in the example scripts used in this package, and deemed to be the
+ is compared to other fits in the example scripts used in this package, and deemed to be the
best reference.
-
+
Args:
T: temperature in kelvin
-
- Reference:
+
+ Reference:
W. Wagner and A. Pruß , "The IAPWS Formulation 1995 for the Thermodynamic Properties
of Ordinary Water Substance for General and Scientific Use", Journal of Physical and Chemical
Reference Data 31, 387-535 (2002) https://doi.org/10.1063/1.1461829
-
+
>>> es_liq(np.asarray([273.16,305.]))
array([ 611.65706974, 4719.32683147])
"""
TvC = constants.temperature_water_vapor_critical_point
PvC = constants.pressure_water_vapor_critical_point
- vt = 1.-T/TvC
- es = PvC * np.exp(TvC/T * (-7.85951783*vt + 1.84408259*vt**1.5 - 11.7866497*vt**3
- + 22.6807411*vt**3.5 - 15.9618719*vt**4 + 1.80122502*vt**7.5))
+ vt = 1.0 - T / TvC
+ es = PvC * np.exp(
+ TvC
+ / T
+ * (
+ -7.85951783 * vt
+ + 1.84408259 * vt**1.5
+ - 11.7866497 * vt**3
+ + 22.6807411 * vt**3.5
+ - 15.9618719 * vt**4
+ + 1.80122502 * vt**7.5
+ )
+ )
return es
+
def es_ice(T):
"""Returns sublimation vapor pressure (Pa) over simple (Ih) ice
-
+
Encodes the emperical fits of Wagner et al., (2011) which also define the IAPWS standard for
sublimation vapor pressure over ice-Ih
-
+
Args:
T: temperature in kelvin
-
+
Reference:
Wagner, W., Riethmann, T., Feistel, R. & Harvey, A. H. New Equations for the Sublimation
Pressure and Melting Pressure of H 2 O Ice Ih. Journal of Physical and Chemical Reference
@@ -78,42 +91,44 @@ def es_ice(T):
TvT = constants.temperature_water_vapor_triple_point
PvT = constants.pressure_water_vapor_triple_point
- a1 = -0.212144006e+2
- a2 = 0.273203819e+2
- a3 = -0.610598130e+1
- b1 = 0.333333333e-2
- b2 = 0.120666667e+1
- b3 = 0.170333333e+1
- theta = T/TvT
- es = PvT * np.exp((a1*theta**b1 + a2 * theta**b2 + a3 * theta**b3)/theta)
+ a1 = -0.212144006e2
+ a2 = 0.273203819e2
+ a3 = -0.610598130e1
+ b1 = 0.333333333e-2
+ b2 = 0.120666667e1
+ b3 = 0.170333333e1
+ theta = T / TvT
+ es = PvT * np.exp((a1 * theta**b1 + a2 * theta**b2 + a3 * theta**b3) / theta)
return es
+
def es_mxd(T):
"""Returns the minimum of the sublimation and saturation vapor pressure
-
+
Calculates both the sublimation vapor pressure over ice Ih using es_ice and that over planar
water using es_liq, and returns the minimum of the two quantities.
-
+
Args:
T: temperature in kelvin
-
+
Returns:
value of es_ice(T) for T < 273.15 and es_liq(T) otherwise
-
+
>>> es_mxd(np.asarray([305.,260.]))
array([4719.32683147, 195.80103377])
"""
- return np.minimum(es_liq(T),es_ice(T))
+ return np.minimum(es_liq(T), es_ice(T))
+
def es_liq_murphykoop(T):
"""Returns saturation vapor pressure (Pa) over liquid water
-
+
Encodes the empirical fit (Eq. 10) of Murphy and Koop (2011) which improves on the Wagner and
Pruß fits for supercooled conditions.
-
+
Args:
T: temperature in kelvin
-
+
Reference:
Murphy, D. M. & Koop, T. Review of the vapour pressures of ice and supercooled water for
atmospheric applications. Q. J. R. Meteorol. Soc. 131, 1539–1565 (2005).
@@ -122,18 +137,21 @@ def es_liq_murphykoop(T):
array([6.11657044e+02, 9.39696372e-07])
"""
- X = np.tanh(0.0415*(T - 218.8)) * (53.878 - 1331.22/T - 9.44523 * np.log(T) + 0.014025*T)
- return np.exp(54.842763 - 6763.22/T - 4.210*np.log(T) + 0.000367*T + X)
+ X = np.tanh(0.0415 * (T - 218.8)) * (
+ 53.878 - 1331.22 / T - 9.44523 * np.log(T) + 0.014025 * T
+ )
+ return np.exp(54.842763 - 6763.22 / T - 4.210 * np.log(T) + 0.000367 * T + X)
+
def es_liq_hardy(T):
"""Returns satruation vapor pressure (Pa) over liquid water
-
- Encodes the empirical fit (Eq. 10) of Hardy (1998) which is often used in the postprocessing
+
+ Encodes the empirical fit (Eq. 10) of Hardy (1998) which is often used in the postprocessing
of radiosondes
-
+
Args:
T: temperature in kelvin
-
+
Reference:
Hardy, B., 1998, ITS-90 Formulations for Vapor Pressure, Frostpoint Temperature, Dewpoint
Temperature, and Enhancement Factors in the Range –100 to +100 °C, The Proceedings of the
@@ -142,61 +160,71 @@ def es_liq_hardy(T):
>>> es_liq_hardy(np.asarray([273.16,260.]))
array([611.65715494, 222.65143353])
"""
- X = (-2.8365744e+3/(T*T) - 6.028076559e+3/T + 19.54263612 - 2.737830188e-2*T
- + 1.6261698e-5*T**2 + 7.0229056e-10*T**3 - 1.8680009e-13*T**4 + 2.7150305 * np.log(T))
+ X = (
+ -2.8365744e3 / (T * T)
+ - 6.028076559e3 / T
+ + 19.54263612
+ - 2.737830188e-2 * T
+ + 1.6261698e-5 * T**2
+ + 7.0229056e-10 * T**3
+ - 1.8680009e-13 * T**4
+ + 2.7150305 * np.log(T)
+ )
return np.exp(X)
+
def es_liq_analytic(T, delta_cl=constants.delta_cl):
"""Analytic approximation for saturation vapor pressure over iquid
-
- Uses the rankine (constant specific heat, negligible condensate volume) approximations to
+
+ Uses the rankine (constant specific heat, negligible condensate volume) approximations to
calculate the saturation vapor pressure over liquid. The procedure is described in Eq(4) of
Romps (2017) and best approximates the actual value for specific heats that differ slightly
- from the best estimates of these quantities which are provided as default quantities.
+ from the best estimates of these quantities which are provided as default quantities.
Romps recommends cl = 4119 J/kg/K, and cpv = 1861 J/kg/K.
-
+
Args:
T: temperature in kelvin
delta_cl: differnce between isobaric specific heat capacity of vapor and that of liquid.
-
+
Returns:
value of saturation vapor pressure over liquid water in Pa
-
+
Reference:
Romps, D. M. Exact Expression for the Lifting Condensation Level. Journal of the Atmospheric
Sciences 74, 3891–3900 (2017).
Romps, D. M. Accurate expressions for the dew point and frost point derived from the Rankine-
Kirchhoff approximations. Journal of the Atmospheric Sciences (2021) doi:10.1175/JAS-D-20-0301.1.
-
+
>>> es_liq_analytic(np.asarray([273.16,305.]))
array([ 611.655 , 4711.13161169])
"""
TvT = constants.temperature_water_vapor_triple_point
PvT = constants.pressure_water_vapor_triple_point
lvT = constants.vaporization_enthalpy_triple_point
- Rv = constants.water_vapor_gas_constant
+ Rv = constants.water_vapor_gas_constant
- c1 = delta_cl/Rv
- c2 = lvT/(Rv*TvT) - c1
- es = PvT * np.exp(c2*(1.-TvT/T)) * (T/TvT)**c1
+ c1 = delta_cl / Rv
+ c2 = lvT / (Rv * TvT) - c1
+ es = PvT * np.exp(c2 * (1.0 - TvT / T)) * (T / TvT) ** c1
return es
+
def es_ice_analytic(T, delta_ci=constants.delta_ci):
"""Analytic approximation for saturation vapor pressure over ice
-
- Uses the rankine (constant specific heat, negligible condensate volume) approximations to
+
+ Uses the rankine (constant specific heat, negligible condensate volume) approximations to
calculate the saturation vapor pressure over ice. The procedure is described in Eq(4) of
Romps (2017) and best approximates the actual value for specific heats that differ slightly
- from the best estimates of these quantities which are provided as default quantities.
+ from the best estimates of these quantities which are provided as default quantities.
Romps recommends ci = 1861 J/kg/K, and cpv = 1879 J/kg/K.
-
+
Args:
T: temperature in kelvin
delta_cl: differnce between isobaric specific heat capacity of vapor and that of liquid.
-
+
Returns:
value of saturation vapor pressure over liquid water in Pa
-
+
Reference:
Romps, D. M. Exact Expression for the Lifting Condensation Level. Journal of the Atmospheric
Sciences 74, 3891–3900 (2017).
@@ -210,38 +238,40 @@ def es_ice_analytic(T, delta_ci=constants.delta_ci):
TvT = constants.temperature_water_vapor_triple_point
PvT = constants.pressure_water_vapor_triple_point
lsT = constants.sublimation_enthalpy_triple_point
- Rv = constants.water_vapor_gas_constant
+ Rv = constants.water_vapor_gas_constant
- c1 = delta_ci/Rv
- c2 = lsT/(Rv*TvT) - c1
- es = PvT * np.exp(c2*(1.-TvT/T)) * (T/TvT)**c1
+ c1 = delta_ci / Rv
+ c2 = lsT / (Rv * TvT) - c1
+ es = PvT * np.exp(c2 * (1.0 - TvT / T)) * (T / TvT) ** c1
return es
+
def es_mxd_analytic(T, delta_cl=constants.delta_cl, delta_ci=constants.delta_ci):
"""Returns the minimum of the analytic sublimation and saturation vapor pressure
-
+
Calculates both the sublimation vapor pressure over ice Ih using es_ice_analytic and
that over planar water using es_liq_analytic, and returns the minimum of the two
quantities.
-
+
Args:
T: temperature in kelvin
-
+
Returns:
value of es_ice_analytic(T) for T < 273.15 and es_liq_analytic(T) otherwise
-
+
>>> es_ice_analytic(np.asarray([273.16,260.]))
array([611.655 , 195.99959431])
"""
- return np.minimum(es_liq_analytic(T,delta_cl),es_ice_analytic(T,delta_ci))
+ return np.minimum(es_liq_analytic(T, delta_cl), es_ice_analytic(T, delta_ci))
-def vaporization_enthalpy(TK,delta_cl=constants.delta_cl):
+
+def vaporization_enthalpy(TK, delta_cl=constants.delta_cl):
"""Returns the vaporization enthlapy of water (J/kg)
-
- The vaporization enthalpy is calculated from a linear depdence on temperature about a
- reference value valid at the melting temperature. This approximation is consistent with the
+
+ The vaporization enthalpy is calculated from a linear depdence on temperature about a
+ reference value valid at the melting temperature. This approximation is consistent with the
assumption of a Rankine fluid.
-
+
Args:
T: temperature in kelvin
delta_cl: differnce between isobaric specific heat capacity of vapor and that of liquid.
@@ -249,17 +279,18 @@ def vaporization_enthalpy(TK,delta_cl=constants.delta_cl):
>>> vaporization_enthalpy(np.asarray([305.,273.15]))
array([2427211.264, 2500930. ])
"""
- T0 = constants.standard_temperature
+ T0 = constants.standard_temperature
lv0 = constants.vaporization_enthalpy_stp
- return lv0 + delta_cl*(TK-T0)
+ return lv0 + delta_cl * (TK - T0)
-def sublimation_enthalpy(TK,delta_ci=constants.delta_ci):
+
+def sublimation_enthalpy(TK, delta_ci=constants.delta_ci):
"""Returns the sublimation enthlapy of water (J/kg)
-
- The sublimation enthalpy is calculated from a linear depdence on temperature about a
- reference value valid at the melting temperature. This approximation is consistent with the
+
+ The sublimation enthalpy is calculated from a linear depdence on temperature about a
+ reference value valid at the melting temperature. This approximation is consistent with the
assumption of a Rankine fluid.
-
+
Args:
T: temperature in kelvin
delta_cl: differnce between isobaric specific heat capacity of vapor and that of liquid.
@@ -268,67 +299,72 @@ def sublimation_enthalpy(TK,delta_ci=constants.delta_ci):
>>> sublimation_enthalpy(273.15)
2834350.0
"""
- T0 = constants.standard_temperature
- ls0 = constants.sublimation_enthalpy_stp
- return ls0 + delta_ci*(TK-T0)
+ T0 = constants.standard_temperature
+ ls0 = constants.sublimation_enthalpy_stp
+ return ls0 + delta_ci * (TK - T0)
-def partial_pressure_to_mixing_ratio(pp,p):
+
+def partial_pressure_to_mixing_ratio(pp, p):
"""Returns the mass mixing ratio given the partial pressure and pressure
-
+
>>> partial_pressure_to_mixing_ratio(es_liq(300.),60000.)
0.0389569254590098
"""
eps1 = constants.rd_over_rv
- return eps1*pp/(p-pp)
+ return eps1 * pp / (p - pp)
+
-def mixing_ratio_to_partial_pressure(r,p):
+def mixing_ratio_to_partial_pressure(r, p):
"""Returns the partial pressure (pp in units of p) from a gas' mixing ratio
-
+
Args:
r: mass mixing ratio (unitless)
- p: pressure in same units as desired return value
-
+ p: pressure in same units as desired return value
+
>>> mixing_ratio_to_partial_pressure(2e-5,60000.)
1.929375975915276
"""
eps1 = constants.rd_over_rv
- return r*p/(eps1+r)
+ return r * p / (eps1 + r)
+
+
+def partial_pressure_to_specific_humidity(pp, p):
+ """Returns the specific mass given the partial pressure and pressure.
-def partial_pressure_to_specific_humidity(pp,p):
- """Returns the specific mass given the partial pressure and pressure.
-
The specific mass can be written in terms of partial pressure and pressure as
- expressed here only if the gas quanta contains no condensate phases. In this
+ expressed here only if the gas quanta contains no condensate phases. In this
case the specific humidity is the same as the co-dryair specific mass. In
- situations where condensate is present one should instead calculate
+ situations where condensate is present one should instead calculate
$q = r*(1-qt)$ which would require an additional argument
>>> partial_pressure_to_specific_humidity(es_liq(300.),60000.)
0.037496189210922945
- """
- r = partial_pressure_to_mixing_ratio(pp,p)
- return r/(1+r)
+ """
+ r = partial_pressure_to_mixing_ratio(pp, p)
+ return r / (1 + r)
-def saturation_partition(P,ps,qt):
+
+def saturation_partition(P, ps, qt):
"""Returns the water vapor specific humidity given saturation vapor presure
-
+
When condensate is present the saturation specific humidity and the total
specific humidity differ, and the latter weights the mixing ratio when
calculating the former from the saturation mixing ratio. In subsaturated air
the vapor speecific humidity is just the total specific humidity
- """
- qs = partial_pressure_to_mixing_ratio(ps,P) * (1. - qt)
- return np.minimum(qt,qs)
-
-def theta(TK,PPa,qv=0., ql=0., qi=0.):
+ """
+ qs = partial_pressure_to_mixing_ratio(ps, P) * (1.0 - qt)
+ return np.minimum(qt, qs)
+
+
+def theta(TK, PPa, qv=0.0, ql=0.0, qi=0.0):
"""Returns the potential temperature for an unsaturated moist fluid
-
+
This expressed the potential temperature in away that makes it possible to account
for the influence of the specific water mass (in different phases) to influence the
adiabatic factor R/cp. The default is the usualy dry potential temperature.
-
+
Args:
TK: temperature in kelvin
PPa: pressure in pascal
@@ -337,450 +373,497 @@ def theta(TK,PPa,qv=0., ql=0., qi=0.):
qi: specific ice mass
es: form of the saturation vapor pressure to use
-
+
"""
- Rd = constants.dry_air_gas_constant
- Rv = constants.water_vapor_gas_constant
- cpd = constants.isobaric_dry_air_specific_heat
- cpv = constants.isobaric_water_vapor_specific_heat
- cl = constants.liquid_water_specific_heat
- ci = constants.frozen_water_specific_heat
- P0 = constants.P0
-
- qd = 1.0-qv-ql-qi
- kappa = (qd*Rd + qv*Rv)/(qd*cpd + qv*cpv + ql*cl + qi*ci)
- return TK*(P0/PPa)**kappa
-
-def theta_e_bolton(TK,PPa,qt,es=es_liq):
+ Rd = constants.dry_air_gas_constant
+ Rv = constants.water_vapor_gas_constant
+ cpd = constants.isobaric_dry_air_specific_heat
+ cpv = constants.isobaric_water_vapor_specific_heat
+ cl = constants.liquid_water_specific_heat
+ ci = constants.frozen_water_specific_heat
+ P0 = constants.P0
+
+ qd = 1.0 - qv - ql - qi
+ kappa = (qd * Rd + qv * Rv) / (qd * cpd + qv * cpv + ql * cl + qi * ci)
+ return TK * (P0 / PPa) ** kappa
+
+
+def theta_e_bolton(TK, PPa, qt, es=es_liq):
"""Returns the pseudo equivalent potential temperature.
-
+
Following Eq. 43 in Bolton (1980) the (pseudo) equivalent potential temperature
is calculated and returned by this function
-
+
Args:
TK: temperature in kelvin
PPa: pressure in pascal
qt: specific total water mass
es: form of the saturation vapor pressure to use
-
+
Reference:
- Bolton, D. The Computation of Equivalent Potential Temperature. Monthly Weather
+ Bolton, D. The Computation of Equivalent Potential Temperature. Monthly Weather
Review 108, 1046–1053 (1980).
"""
- P0 = constants.standard_pressure
- p2r = partial_pressure_to_mixing_ratio
- r2p = mixing_ratio_to_partial_pressure
-
- rv = np.minimum(qt/(1.-qt), p2r(es(TK),PPa)) # mixing ratio of vapor (not gas Rv)
- pv = r2p(rv,PPa)
+ P0 = constants.standard_pressure
+ p2r = partial_pressure_to_mixing_ratio
+ r2p = mixing_ratio_to_partial_pressure
- TL = 55.0 + 2840./(3.5*np.log(TK) - np.log(pv/100.) - 4.805)
- return TK*(P0/PPa)**(0.2854*(1.0 - 0.28*rv)) * np.exp((3376./TL - 2.54)*rv*(1+0.81*rv))
+ rv = np.minimum(
+ qt / (1.0 - qt), p2r(es(TK), PPa)
+ ) # mixing ratio of vapor (not gas Rv)
+ pv = r2p(rv, PPa)
-def theta_e(TK,PPa,qt,es=es_liq):
+ TL = 55.0 + 2840.0 / (3.5 * np.log(TK) - np.log(pv / 100.0) - 4.805)
+ return (
+ TK
+ * (P0 / PPa) ** (0.2854 * (1.0 - 0.28 * rv))
+ * np.exp((3376.0 / TL - 2.54) * rv * (1 + 0.81 * rv))
+ )
+
+
+def theta_e(TK, PPa, qt, es=es_liq):
"""Returns the equivalent potential temperature
-
+
Follows Eq. 11 in Marquet and Stevens (2022). The closed form solutionis derived for a
Rankine-Kirchoff fluid (constant specific heats). Differences arising from its
- calculation using more accurate expressions (such as the default) as opposed to less
+ calculation using more accurate expressions (such as the default) as opposed to less
accurate, but more consistent, formulations are on the order of millikelvin
-
+
Args:
TK: temperature in kelvin
PPa: pressure in pascal
qt: total water specific humidity (unitless)
es: form of the saturation vapor pressure
-
+
Reference:
Marquet, P. & Stevens, B. On Moist Potential Temperatures and Their Ability to
Characterize Differences in the Properties of Air Parcels. Journal of the Atmospheric
Sciences 79, 1089–1103 (2022).
"""
- P0 = constants.standard_pressure
- Rd = constants.dry_air_gas_constant
- Rv = constants.water_vapor_gas_constant
- cpd = constants.isobaric_dry_air_specific_heat
- cl = constants.liquid_water_specific_heat
- lv = vaporization_enthalpy
+ P0 = constants.standard_pressure
+ Rd = constants.dry_air_gas_constant
+ Rv = constants.water_vapor_gas_constant
+ cpd = constants.isobaric_dry_air_specific_heat
+ cl = constants.liquid_water_specific_heat
+ lv = vaporization_enthalpy
ps = es(TK)
- qv = saturation_partition(PPa,ps,qt)
-
- Re = (1.0-qt)*Rd
- R = Re + qv*Rv
- pv = qv * (Rv/R) *PPa
- RH = pv/ps
- cpe= cpd + qt*(cl-cpd)
- omega_e = RH**(-qv*Rv/cpe) * (R/Re)**(Re/cpe)
- theta_e = TK*(P0/PPa)**(Re/cpe)*omega_e*np.exp(qv*lv(TK)/(cpe*TK))
- return(theta_e)
-
-def theta_l(TK,PPa,qt,es=es_liq):
+ qv = saturation_partition(PPa, ps, qt)
+
+ Re = (1.0 - qt) * Rd
+ R = Re + qv * Rv
+ pv = qv * (Rv / R) * PPa
+ RH = pv / ps
+ cpe = cpd + qt * (cl - cpd)
+ omega_e = RH ** (-qv * Rv / cpe) * (R / Re) ** (Re / cpe)
+ theta_e = TK * (P0 / PPa) ** (Re / cpe) * omega_e * np.exp(qv * lv(TK) / (cpe * TK))
+ return theta_e
+
+
+def theta_l(TK, PPa, qt, es=es_liq):
"""Returns the liquid-water potential temperature
-
+
Follows Eq. 16 in Marquet and Stevens (2022). The closed form solutionis derived for a
Rankine-Kirchoff fluid (constant specific heats). Differences arising from its
- calculation using more accurate expressions (such as the default) as opposed to less
+ calculation using more accurate expressions (such as the default) as opposed to less
accurate, but more consistent, formulations are on the order of millikelvin
-
+
Args:
TK: temperature in kelvin
PPa: pressure in pascal
qt: total water specific humidity (unitless)
es: form of the saturation vapor pressure
-
+
Reference:
Marquet, P. & Stevens, B. On Moist Potential Temperatures and Their Ability to
Characterize Differences in the Properties of Air Parcels. Journal of the Atmospheric
Sciences 79, 1089–1103 (2022).
"""
- P0 = constants.standard_pressure
- Rd = constants.dry_air_gas_constant
- Rv = constants.water_vapor_gas_constant
- cpd = constants.isobaric_dry_air_specific_heat
- cpv = constants.isobaric_water_vapor_specific_heat
- lv = vaporization_enthalpy
+ P0 = constants.standard_pressure
+ Rd = constants.dry_air_gas_constant
+ Rv = constants.water_vapor_gas_constant
+ cpd = constants.isobaric_dry_air_specific_heat
+ cpv = constants.isobaric_water_vapor_specific_heat
+ lv = vaporization_enthalpy
ps = es(TK)
- qv = saturation_partition(PPa,ps,qt)
- ql = qt-qv
-
- R = Rd*(1-qt) + qv*Rv
- Rl = Rd + qt*(Rv - Rd)
- cpl= cpd + qt*(cpv-cpd)
-
- omega_l = (R/Rl)**(Rl/cpl) * (qt/(qv+1.e-15))**(qt*Rv/cpl)
- theta_l = (TK*(P0/PPa)**(Rl/cpl)) *omega_l*np.exp(-ql*lv(TK)/(cpl*TK))
- return(theta_l)
-
-def theta_s(TK,PPa,qt,es=es_liq):
+ qv = saturation_partition(PPa, ps, qt)
+ ql = qt - qv
+
+ R = Rd * (1 - qt) + qv * Rv
+ Rl = Rd + qt * (Rv - Rd)
+ cpl = cpd + qt * (cpv - cpd)
+
+ omega_l = (R / Rl) ** (Rl / cpl) * (qt / (qv + 1.0e-15)) ** (qt * Rv / cpl)
+ theta_l = (
+ (TK * (P0 / PPa) ** (Rl / cpl)) * omega_l * np.exp(-ql * lv(TK) / (cpl * TK))
+ )
+ return theta_l
+
+
+def theta_s(TK, PPa, qt, es=es_liq):
"""Returns the entropy potential temperature
-
+
Follows Eq. 18 in Marquet and Stevens (2022). The closed form solutionis derived for a
Rankine-Kirchoff fluid (constant specific heats). Differences arising from its
- calculation using more accurate expressions (such as the default) as opposed to less
+ calculation using more accurate expressions (such as the default) as opposed to less
accurate, but more consistent, formulations are on the order of millikelvin
-
+
Args:
TK: temperature in kelvin
PPa: pressure in pascal
qt: total water specific humidity (unitless)
es: form of the saturation vapor pressure
-
+
Reference:
Marquet, P. & Stevens, B. On Moist Potential Temperatures and Their Ability to
Characterize Differences in the Properties of Air Parcels. Journal of the Atmospheric
Sciences 79, 1089–1103 (2022).
-
+
Marquet, P. Definition of a moist entropy potential temperature: application to FIRE-I
data flights: Moist Entropy Potential Temperature. Q.J.R. Meteorol. Soc. 137, 768–791 (2011).
"""
- P0 = constants.standard_pressure
- T0 = constants.standard_temperature
+ P0 = constants.standard_pressure
+ T0 = constants.standard_temperature
sd00 = constants.entropy_dry_air_satmt
sv00 = constants.entropy_water_vapor_satmt
- Rd = constants.dry_air_gas_constant
- Rv = constants.water_vapor_gas_constant
- cpd = constants.isobaric_dry_air_specific_heat
- cpv = constants.isobaric_water_vapor_specific_heat
+ Rd = constants.dry_air_gas_constant
+ Rv = constants.water_vapor_gas_constant
+ cpd = constants.isobaric_dry_air_specific_heat
+ cpv = constants.isobaric_water_vapor_specific_heat
eps1 = constants.rd_over_rv
eps2 = constants.rv_over_rd_minus_one
- lv = vaporization_enthalpy
+ lv = vaporization_enthalpy
- kappa = Rd/cpd
- e0 = es(T0)
- Lmbd = ((sv00 - Rv*np.log(e0/P0)) - (sd00 - Rd*np.log(1-e0/P0)))/cpd
- lmbd = cpv/cpd - 1.
- eta = 1/eps1
+ kappa = Rd / cpd
+ e0 = es(T0)
+ Lmbd = ((sv00 - Rv * np.log(e0 / P0)) - (sd00 - Rd * np.log(1 - e0 / P0))) / cpd
+ lmbd = cpv / cpd - 1.0
+ eta = 1 / eps1
delta = eps2
- gamma = kappa/eps1
- r0 = e0/(P0-e0)/eta
+ gamma = kappa / eps1
+ r0 = e0 / (P0 - e0) / eta
ps = es(TK)
- qv = saturation_partition(PPa,ps,qt)
- ql = qt-qv
-
- R = Rd + qv*(Rv - Rd)
- pv = qv * (Rv/R) *PPa
- RH = pv/ps
- rv = qv/(1-qv)
-
- x1 = (TK/T0)**(lmbd*qt) * (P0/PPa)**(kappa*delta*qt) * (rv/r0)**(-gamma*qt) * RH**(gamma*ql)
- x2 = (1.+eta*rv)**(kappa*(1.+delta*qt)) * (1.+eta*r0)**(-kappa*delta*qt)
- theta_s = (TK*(P0/PPa)**(kappa)) * np.exp(-ql*lv(TK)/(cpd*TK)) * np.exp(qt*Lmbd) * x1 * x2
- return(theta_s)
-
-def theta_es(TK,PPa,es=es_liq):
+ qv = saturation_partition(PPa, ps, qt)
+ ql = qt - qv
+
+ R = Rd + qv * (Rv - Rd)
+ pv = qv * (Rv / R) * PPa
+ RH = pv / ps
+ rv = qv / (1 - qv)
+
+ x1 = (
+ (TK / T0) ** (lmbd * qt)
+ * (P0 / PPa) ** (kappa * delta * qt)
+ * (rv / r0) ** (-gamma * qt)
+ * RH ** (gamma * ql)
+ )
+ x2 = (1.0 + eta * rv) ** (kappa * (1.0 + delta * qt)) * (1.0 + eta * r0) ** (
+ -kappa * delta * qt
+ )
+ theta_s = (
+ (TK * (P0 / PPa) ** (kappa))
+ * np.exp(-ql * lv(TK) / (cpd * TK))
+ * np.exp(qt * Lmbd)
+ * x1
+ * x2
+ )
+ return theta_s
+
+
+def theta_es(TK, PPa, es=es_liq):
"""Returns the saturated equivalent potential temperature
-
+
Adapted from Eq. 11 in Marquet and Stevens (2022) with the assumption that the gas quanta is
everywhere just saturated.
-
+
Args:
TK: temperature in kelvin
PPa: pressure in pascal
qt: total water specific humidity (unitless)
es: form of the saturation vapor pressure
-
+
Reference:
Characterize Differences in the Properties of Air Parcels. Journal of the Atmospheric
Sciences 79, 1089–1103 (2022).
"""
- P0 = constants.standard_pressure
- Rd = constants.dry_air_gas_constant
- Rv = constants.water_vapor_gas_constant
- cpd = constants.isobaric_dry_air_specific_heat
- cl = constants.liquid_water_specific_heat
- p2q = partial_pressure_to_specific_humidity
- lv = vaporization_enthalpy
+ P0 = constants.standard_pressure
+ Rd = constants.dry_air_gas_constant
+ Rv = constants.water_vapor_gas_constant
+ cpd = constants.isobaric_dry_air_specific_heat
+ cl = constants.liquid_water_specific_heat
+ p2q = partial_pressure_to_specific_humidity
+ lv = vaporization_enthalpy
ps = es(TK)
- qs = p2q(ps,PPa)
+ qs = p2q(ps, PPa)
+
+ Re = (1.0 - qs) * Rd
+ R = Re + qs * Rv
+ cpe = cpd + qs * (cl - cpd)
+ omega_e = (R / Re) ** (Re / cpe)
+ theta_es = (
+ TK * (P0 / PPa) ** (Re / cpe) * omega_e * np.exp(qs * lv(TK) / (cpe * TK))
+ )
+ return theta_es
- Re = (1.0-qs)*Rd
- R = Re + qs*Rv
- cpe= cpd + qs*(cl-cpd)
- omega_e = (R/Re)**(Re/cpe)
- theta_es = TK*(P0/PPa)**(Re/cpe)*omega_e*np.exp(qs*lv(TK)/(cpe*TK))
- return(theta_es)
-def theta_rho(TK,PPa,qt,es=es_liq):
+def theta_rho(TK, PPa, qt, es=es_liq):
"""Returns the density liquid-water potential temperature
-
+
calculates $\theta_\mathrm{l} R/R_\mathrm{d}$ where $R$ is the gas constant of a
most fluid. For an unsaturated fluid this is identical to the density potential
temperature baswed on the two component fluid thermodynamic constants.
-
+
Args:
TK: temperature in kelvin
PPa: pressure in pascal
qt: total water specific humidity (unitless)
es: form of the saturation vapor pressure
"""
- Rd = constants.dry_air_gas_constant
- Rv = constants.water_vapor_gas_constant
+ Rd = constants.dry_air_gas_constant
+ Rv = constants.water_vapor_gas_constant
ps = es(TK)
- qv = saturation_partition(PPa,ps,qt)
- theta_rho = theta_l(TK,PPa,qt,es) * (1.-qt + qv*Rv/Rd)
- return(theta_rho)
+ qv = saturation_partition(PPa, ps, qt)
+ theta_rho = theta_l(TK, PPa, qt, es) * (1.0 - qt + qv * Rv / Rd)
+ return theta_rho
+
-def T_from_Te(Te,P,qt,es=es_liq):
+def T_from_Te(Te, P, qt, es=es_liq):
"""Returns temperature for an atmosphere whose state is given by theta_e
-
- This equation allows the temperature to be inferred from a state description
- in terms of theta_e. It derives temperature by numerically inverting the
- expression for theta_e provided in this package. Uses a least-squares non-linear
- optimization to find the value of T such that $theta_e(T,P,q) = theta_e$
-
-Args:
- Te: equivalent potential temperature in kelvin
- P: pressure in pascal
- qt: total water specific humidity (unitless)
- es: form of the saturation vapor pressure
- >>> T_from_Te(350.,100000.,17.e-3)
- array([304.49321301])
+ This equation allows the temperature to be inferred from a state description
+ in terms of theta_e. It derives temperature by numerically inverting the
+ expression for theta_e provided in this package. Uses a least-squares non-linear
+ optimization to find the value of T such that $theta_e(T,P,q) = theta_e$
+
+ Args:
+ Te: equivalent potential temperature in kelvin
+ P: pressure in pascal
+ qt: total water specific humidity (unitless)
+ es: form of the saturation vapor pressure
+
+ >>> T_from_Te(350.,100000.,17.e-3)
+ array([304.49321301])
"""
- def zero(T,Te,P,qt,es):
- return np.abs(Te-theta_e(T,P,qt,es))
- return optimize.fsolve(zero, 280., args=(Te,P,qt,es), xtol=1.e-10)
-def T_from_Tl(Tl,P,qt,es=es_liq):
+ def zero(T, Te, P, qt, es):
+ return np.abs(Te - theta_e(T, P, qt, es))
+
+ return optimize.fsolve(zero, 280.0, args=(Te, P, qt, es), xtol=1.0e-10)
+
+
+def T_from_Tl(Tl, P, qt, es=es_liq):
"""returns temperature for an atmosphere whose state is given by theta_l
-
- This equation allows the temperature to be inferred from a state description
- in terms of theta_l. It derives temperature by numerically inverting the
- expression for theta_l provided in this package. Uses a least-squares non-linear
- optimization to find the value of T such that $theta_l(T,P,q) = theta_l$
-Args:
- Tl: liquid-water potential temperature in kelvin
- P: pressure in pascal
- qt: total water specific humidity (unitless)
- es: form of the saturation vapor pressure
+ This equation allows the temperature to be inferred from a state description
+ in terms of theta_l. It derives temperature by numerically inverting the
+ expression for theta_l provided in this package. Uses a least-squares non-linear
+ optimization to find the value of T such that $theta_l(T,P,q) = theta_l$
+
+ Args:
+ Tl: liquid-water potential temperature in kelvin
+ P: pressure in pascal
+ qt: total water specific humidity (unitless)
+ es: form of the saturation vapor pressure
- >>> T_from_Tl(282., 90000, 20.e-3)
- array([289.73684039])
+ >>> T_from_Tl(282., 90000, 20.e-3)
+ array([289.73684039])
"""
- def zero(T,Tl,P,qt,es):
- return np.abs(Tl-theta_l(T,P,qt,es))
-
- return optimize.fsolve(zero, 280., args=(Tl,P,qt,es), xtol=1.e-10)
-def T_from_Ts(Ts,P,qt,es=es_liq):
+ def zero(T, Tl, P, qt, es):
+ return np.abs(Tl - theta_l(T, P, qt, es))
+
+ return optimize.fsolve(zero, 280.0, args=(Tl, P, qt, es), xtol=1.0e-10)
+
+
+def T_from_Ts(Ts, P, qt, es=es_liq):
"""Returns temperature for an atmosphere whose state is given by theta_s
-
- This equation allows the temperature to be inferred from a state description
- in terms of theta_s. It derives temperature by numerically inverting the
- expression for theta_s provided in this package. Uses a least-squares non-linear
- optimization to find the value of T such that $theta_s(T,P,q) = theta_s$
-
-Args:
- Ts: entropy potential temperature in kelvin
- P: pressure in pascal
- qt: total water specific humidity (unitless)
- es: form of the saturation vapor pressure
- >>> T_from_Tl(282.75436951,90000,20.e-3)
- array([289.98864293])
+ This equation allows the temperature to be inferred from a state description
+ in terms of theta_s. It derives temperature by numerically inverting the
+ expression for theta_s provided in this package. Uses a least-squares non-linear
+ optimization to find the value of T such that $theta_s(T,P,q) = theta_s$
+
+ Args:
+ Ts: entropy potential temperature in kelvin
+ P: pressure in pascal
+ qt: total water specific humidity (unitless)
+ es: form of the saturation vapor pressure
+
+ >>> T_from_Tl(282.75436951,90000,20.e-3)
+ array([289.98864293])
"""
- def zero(T,Ts,P,qt,es):
- return np.abs(Ts-theta_s(T,P,qt,es))
-
- return optimize.fsolve(zero, 280., args=(Ts,P,qt,es), xtol=1.e-10)
-def P_from_Te(Te,T,qt,es=es_liq):
+ def zero(T, Ts, P, qt, es):
+ return np.abs(Ts - theta_s(T, P, qt, es))
+
+ return optimize.fsolve(zero, 280.0, args=(Ts, P, qt, es), xtol=1.0e-10)
+
+
+def P_from_Te(Te, T, qt, es=es_liq):
"""Returns pressure for an atmosphere whose state is given by theta_e
-
- This equation allows the pressure to be inferred from a state description
- in terms of theta_e. It derives pressure by numerically inverting the
- expression for theta_e provided in this package. Uses a least-squares non-linear
- optimization to find the value of P such that $theta_e(T,P,q) = theta_e$
-
-Args:
- Tl: liquid-water potential temperature in kelvin
- T: temperature in kelvin
- qt: total water specific humidity (unitless)
- es: form of the saturation vapor pressure
- >>> P_from_Te(350.,305.,17e-3)
- array([100586.3357635])
+ This equation allows the pressure to be inferred from a state description
+ in terms of theta_e. It derives pressure by numerically inverting the
+ expression for theta_e provided in this package. Uses a least-squares non-linear
+ optimization to find the value of P such that $theta_e(T,P,q) = theta_e$
+
+ Args:
+ Tl: liquid-water potential temperature in kelvin
+ T: temperature in kelvin
+ qt: total water specific humidity (unitless)
+ es: form of the saturation vapor pressure
+
+ >>> P_from_Te(350.,305.,17e-3)
+ array([100586.3357635])
"""
- def zero(P,Te,T,qt,es):
- return np.abs(Te-theta_e(T,P,qt,es))
-
- return optimize.fsolve(zero, 90000., args=(Te,T,qt,es), xtol=1.e-10)
-def P_from_Tl(Tl,T,qt,es=es_liq):
+ def zero(P, Te, T, qt, es):
+ return np.abs(Te - theta_e(T, P, qt, es))
+
+ return optimize.fsolve(zero, 90000.0, args=(Te, T, qt, es), xtol=1.0e-10)
+
+
+def P_from_Tl(Tl, T, qt, es=es_liq):
"""Returns pressure for an atmosphere whose state is given by theta_l
-
+
This equation allows the pressure to be inferred from a state description
in terms of theta_l. It derives pressure by numerically inverting the
expression for theta_l provided in this package. Uses a least-squares non-linear
optimization to find the value of P such that $theta_l(T,P,q) = theta_l$
-
+
Args:
Tl: liquid-water potential temperature in kelvin
T: temperature in kelvin
qt: total water specific humidity (unitless)
es: form of the saturation vapor pressure
- >>> P_from_Tl(282.75436951,290,20.e-3)
- array([90027.65146427])
+ >>> P_from_Tl(282.75436951,290,20.e-3)
+ array([90027.65146427])
"""
- def zero(P,Tl,T,qt,es):
- return np.abs(Tl-theta_l(T,P,qt,es))
-
- return optimize.fsolve(zero, 90000., args=(Tl,T,qt,es), xtol=1.e-10)
-def plcl(TK,PPa,qt,es=es_liq):
+ def zero(P, Tl, T, qt, es):
+ return np.abs(Tl - theta_l(T, P, qt, es))
+
+ return optimize.fsolve(zero, 90000.0, args=(Tl, T, qt, es), xtol=1.0e-10)
+
+
+def plcl(TK, PPa, qt, es=es_liq):
"""Returns the pressure at the lifting condensation level
-
+
Calculates the lifting condensation level pressure using an interative solution under the
constraint of constant theta-l. Exact to within the accuracy of the expression of theta-l
which depends on the expression for the saturation vapor pressure
-
+
Args:
TK: temperature in kelvin
PPa: pressure in pascal
qt: specific total water mass
-
- >>> plcl(300.,102000.,17e-3)
- array([95971.69750248])
+
+ >>> plcl(300.,102000.,17e-3)
+ array([95971.69750248])
"""
-
- def zero(P,Tl,qt,es):
- p2r = partial_pressure_to_mixing_ratio
- TK = T_from_Tl(Tl,P,qt)
- qs = p2r(es(TK),P) * (1. - qt)
- return np.abs(qs/qt-1.)
- Tl = theta_l(TK,PPa,qt,es)
- return optimize.fsolve(zero, 80000., args=(Tl,qt,es), xtol=1.e-5)
+ def zero(P, Tl, qt, es):
+ p2r = partial_pressure_to_mixing_ratio
+ TK = T_from_Tl(Tl, P, qt)
+ qs = p2r(es(TK), P) * (1.0 - qt)
+ return np.abs(qs / qt - 1.0)
+
+ Tl = theta_l(TK, PPa, qt, es)
+ return optimize.fsolve(zero, 80000.0, args=(Tl, qt, es), xtol=1.0e-5)
-def plcl_bolton(TK,PPa,qt):
+
+def plcl_bolton(TK, PPa, qt):
"""Returns the pressure at the lifting condensation level
-
+
Following Bolton (1980) the lifting condensation level pressure is derived from the state
of an air parcel. Usually accurate to within about 10 Pa, or about 1 m
-
+
Args:
TK: temperature in kelvin
PPa: pressure in pascal
qt: specific total water mass
-
+
Reference:
- Bolton, D. The Computation of Equivalent Potential Temperature. Monthly Weather
+ Bolton, D. The Computation of Equivalent Potential Temperature. Monthly Weather
Review 108, 1046–1053 (1980).
-
- >>> plcl_bolton(300.,102000.,17e-3)
- 95980.41895404423
+
+ >>> plcl_bolton(300.,102000.,17e-3)
+ 95980.41895404423
"""
- Rd = constants.dry_air_gas_constant
- Rv = constants.water_vapor_gas_constant
- cpd = constants.isobaric_dry_air_specific_heat
- cpv = constants.isobaric_water_vapor_specific_heat
- r2p = mixing_ratio_to_partial_pressure
-
- cp = cpd + qt*(cpv-cpd)
- R = Rd + qt*(Rv-Rd)
- pv = r2p(qt/(1.-qt),PPa)
- Tl = 55 + 2840./(3.5*np.log(TK) - np.log(pv/100.) - 4.805)
- return PPa * (Tl/TK)**(cp/R)
-
-def zlcl(Plcl,T,P,qt,z):
+ Rd = constants.dry_air_gas_constant
+ Rv = constants.water_vapor_gas_constant
+ cpd = constants.isobaric_dry_air_specific_heat
+ cpv = constants.isobaric_water_vapor_specific_heat
+ r2p = mixing_ratio_to_partial_pressure
+
+ cp = cpd + qt * (cpv - cpd)
+ R = Rd + qt * (Rv - Rd)
+ pv = r2p(qt / (1.0 - qt), PPa)
+ Tl = 55 + 2840.0 / (3.5 * np.log(TK) - np.log(pv / 100.0) - 4.805)
+ return PPa * (Tl / TK) ** (cp / R)
+
+
+def zlcl(Plcl, T, P, qt, z):
"""Returns the height of the LCL above mean sea-level
-
+
Given the Plcl, calculate its height in meters given the height of the ambient state
from which it (Plcl) was calculated. This is accomplished by assuming temperature
- changes following a dry adiabat with vertical displacements between the ambient
+ changes following a dry adiabat with vertical displacements between the ambient
temperature and the ambient LCL
-
+
Args:
Plcl: lifting condensation level in Pa
T: ambient temperature in kelvin
P: ambient pressure in pascal
qt: specific total water mass
- z: height at ambient temperature and pressure
-
- >>> zlcl(95000.,300.,90000.,17.e-3,500.)
- 16.621174077862747
+ z: height at ambient temperature and pressure
+
+ >>> zlcl(95000.,300.,90000.,17.e-3,500.)
+ 16.621174077862747
"""
- Rd = constants.dry_air_gas_constant
- Rv = constants.water_vapor_gas_constant
- cpd = constants.isobaric_dry_air_specific_heat
- cpv = constants.isobaric_water_vapor_specific_heat
- g = constants.gravity_earth
+ Rd = constants.dry_air_gas_constant
+ Rv = constants.water_vapor_gas_constant
+ cpd = constants.isobaric_dry_air_specific_heat
+ cpv = constants.isobaric_water_vapor_specific_heat
+ g = constants.gravity_earth
+
+ cp = cpd + qt * (cpv - cpd)
+ R = Rd + qt * (Rv - Rd)
+ return T * (1.0 - (Plcl / P) ** (R / cp)) * cp / g + z
- cp = cpd + qt*(cpv-cpd)
- R = Rd + qt*(Rv-Rd)
- return T*(1. - (Plcl/P)**(R/cp)) * cp / g + z
from scipy.integrate import ode
-def moist_adiabat(Tbeg,Pbeg,Pend,dP,qt,cc=constants.cl,l=vaporization_enthalpy,es = es_liq):
+
+def moist_adiabat(
+ Tbeg, Pbeg, Pend, dP, qt, cc=constants.cl, l=vaporization_enthalpy, es=es_liq
+):
"""Returns the temperature and pressure by integrating along a moist adiabat
-
+
Deriving the moist adiabats by assuming a constant moist potential temperature
provides a Rankine-Kirchoff approximation to the moist adiabat. If thermodynamic
constants are allowed to vary with temperature then the intergation must be
- performed numerically, as outlined here for the case of constant thermodynamic
+ performed numerically, as outlined here for the case of constant thermodynamic
constants and no accounting for the emergence of a solid condensage phase (ice).
-
+
The introduction of this function allows one to estimate, for instance, the effect of
isentropic freezing on the moist adiabat as follows:
-
+
Tliq,Px= moist_adiabat(Tsfc,Psfc,Ptop,dP,qt,cc=constants.cl,l=mt.vaporization_enthalpy,es = mt.es_mxd)
- Tice,Py= moist_adiabat(Tsfc,Psfc,Ptop,dP,qt,cc=constants.ci,l=mt.sublimation_enthalpy ,es = mt.es_mxd)
-
+ Tice,Py= moist_adiabat(Tsfc,Psfc,Ptop,dP,qt,cc=constants.ci,l=mt.sublimation_enthalpy ,es = mt.es_mxd)
+
T = np.ones(len(Tx))*constants.T0
T[Tliq>constants.T0] = Tliq[Tliq>constants.T0]
T[Tice<constants.T0] = Tice[Tice<constants.T0]
-
+
which introduces an isothermal layer in the region where the fusion enthalpy is sufficient to do
the expansional work
-
+
Args:
Tbeg: temperature at P0 in kelvin
Pbeg: starting pressure in pascal
@@ -788,39 +871,40 @@ def moist_adiabat(Tbeg,Pbeg,Pend,dP,qt,cc=constants.cl,l=vaporization_enthalpy,e
dP: integration step
qt: specific mass of total water
es: saturation vapor expression
-
+
"""
- def f(P,T,qt,cc,l):
- Rd = constants.Rd
- Rv = constants.Rv
+
+ def f(P, T, qt, cc, l):
+ Rd = constants.Rd
+ Rv = constants.Rv
cpd = constants.cpd
cpv = constants.cpv
-
- qv = saturation_partition(P,es(T),qt)
- qc = qt-qv
- qd = 1.-qt
-
- R = qd * Rd + qv * Rv
- cp = qd * cpd + qv * cpv + qc * cc
- vol = R * T/P
-
- dX_dT = cp
- dX_dP = vol
- if (qc > 0.):
- beta_P = R/(qd*Rd)
- beta_T = beta_P * l(T)/(Rv * T)
-
- dX_dT += l(T) * qv * beta_T/T
- dX_dP *= ( 1.0 + l(T) * qv * beta_P/(R*T))
- return dX_dP/dX_dT;
-
+
+ qv = saturation_partition(P, es(T), qt)
+ qc = qt - qv
+ qd = 1.0 - qt
+
+ R = qd * Rd + qv * Rv
+ cp = qd * cpd + qv * cpv + qc * cc
+ vol = R * T / P
+
+ dX_dT = cp
+ dX_dP = vol
+ if qc > 0.0:
+ beta_P = R / (qd * Rd)
+ beta_T = beta_P * l(T) / (Rv * T)
+
+ dX_dT += l(T) * qv * beta_T / T
+ dX_dP *= 1.0 + l(T) * qv * beta_P / (R * T)
+ return dX_dP / dX_dT
+
Tx = []
Px = []
- r = ode(f).set_integrator('lsoda',atol=0.0001)
- r.set_initial_value(Tbeg, Pbeg).set_f_params(qt,cc,l)
+ r = ode(f).set_integrator("lsoda", atol=0.0001)
+ r.set_initial_value(Tbeg, Pbeg).set_f_params(qt, cc, l)
while r.successful() and r.t > Pend:
- r.integrate(r.t+dP)
+ r.integrate(r.t + dP)
Tx.append(r.y[0])
Px.append(r.t)
- return np.asarray(Tx),np.asarray(Px)
\ No newline at end of file
+ return np.asarray(Tx), np.asarray(Px)