This notebook is meant to work with MOM6 run via the CESM framework <div class="alert alert-info>
This notebook checks if the fluxes of water, heat, and salt crossing the surface of the ocean are consistent with corresponding changes in the ocean interior. This particular analysis is for ag.e30_beta02.ctrl CESM G_JRA-compset without OBCs.
Results from this notebook:
- Maps of surface boundary fluxes of water, heat, and salt crossing into the liquid seawater in MOM6;
- Computation of self-consistency checks, including global heat, salt, and mass budgets to verify that the model is conserving scalars over the global/regional domain.
Hopes for the use of this notebook:
- To provide a starting point to document boundary fluxes of scalar fields;
- To teach MOM6 users about the boundary fluxes, their patterns, units, and sign conventions;
- To perform self-consistency checks to ensure the model is conserving scalar fields;
- To illustrate a self-contained notebook of use for MOM6 analysis.
This notebook was originally developed at NOAA/GFDL, and it is provided freely to the MOM6 community. GFDL scientists developing MOM6 make extensive use of Python for diagnostics. We solicit modifications/fixes that are useful to the MOM6 community.
# %pylab inline
import matplotlib.pyplot as plt
import netCDF4
import numpy path = '/glade/derecho/scratch/gmarques/b.e30_a06c.cmeps547_GRIS-EVOLVE/run/'
static = netCDF4.Dataset(path+'ocean.mom6.static.nc')
forcing = netCDF4.Dataset(path+'ocean.mom6.frc.0001.nc')
surface = netCDF4.Dataset(path+'ocean.mom6.sfc.0001.nc') ; tvar = 'time'
cp = 3992.0
rho0 = 1035.0
n = 3
print('Variables saved:')
for v in forcing.variables: print (v),
for v in surface.variables: print (v),
for v in static.variables: print (v),Variables saved:
xq
yh
time
nbnd
xh
yq
scalar_axis
taux
tauy
ustar
PRCmE
lprec
fprec
evap
vprec
lrunoff
frunoff
lrunoff_glc
frunoff_glc
seaice_melt
net_heat_coupler
net_heat_surface
frazil
sensible
latent
LwLatSens
SW
LW
Heat_PmE
seaice_melt_heat
heat_content_lrunoff
heat_content_frunoff
heat_content_lprec
heat_content_fprec
heat_content_vprec
heat_content_cond
heat_content_evap
heat_content_frunoff_glc
heat_content_lrunoff_glc
heat_content_massout
heat_content_massin
heat_content_surfwater
internal_heat
heat_added
hfds
p_surf
salt_flux
salt_flux_in
vprec_global_adjustment
net_fresh_water_global_adjustment
salt_flux_global_restoring_adjustment
net_massout
net_massin
average_T1
average_T2
average_DT
time_bounds
xh
yh
time
xq
yq
SSH
SST
SSS
SSU
SSV
mass_wt
temp_int
salt_int
KPP_OBLdepth
xh
yh
time
xq
yq
geolon
geolat
geolon_c
geolat_c
geolon_u
geolat_u
geolon_v
geolat_v
area_t
depth_ocean
wet
wet_c
wet_u
wet_v
Coriolis
# This section fills the fields used in this notebook.
#--------------------------------------------------------------
# geometric factors
lon = static.variables['geolon'][:]
lat = static.variables['geolat'][:]
wet = static.variables['wet'][:]
area = static.variables['area_t'][:]*wet
#--------------------------------------------------------------
# time in days
time = surface.variables[tvar][:]*86400.
#--------------------------------------------------------------
# sea surface temperature
sst = surface.variables['SST'][n]
#--------------------------------------------------------------
# \int \rho \, dz \, (1,Theta,S), with \rho = \rho_{0} for Bousssinesq.
mass_wt = surface.variables['mass_wt']
tomint = surface.variables['temp_int']
somint = surface.variables['salt_int']
#--------------------------------------------------------------
# mass flux of water crossing ocean surface [kg/(m^2 s)]
# positive values indicate mass entering ocean;
# negative values indicate mass leaving ocean.
# net mass flux entering ocean
net_massin = forcing.variables['net_massin'][n]
# net mass flux leaving ocean
net_massout = forcing.variables['net_massout'][n]
# evaporation (negative) and condensation (positive)
evap = forcing.variables['evap'][n]
# meltw, fresh water from melting/forming sea ice
meltw = forcing.variables['seaice_melt'][n]
# liquid runoff entering ocean (non-negative)
lrunoff = forcing.variables['lrunoff'][n]
# frozen runoff entering ocean (non-negative)
frunoff = forcing.variables['frunoff'][n]
# liquid runoff entering ocean (non-negative)
lrunoff_glc = forcing.variables['lrunoff_glc'][n]
# frozen runoff entering ocean (non-negative)
frunoff_glc = forcing.variables['frunoff_glc'][n]
# liquid precipitation entering ocean.
# note: includes exchanges with sea-ice, with
# melt adding mass to ocean; formation removing mass.
lprec = forcing.variables['lprec'][n]
# frozen precipitation entering ocean.
fprec = forcing.variables['fprec'][n]
# virtual precipitation arising from conversion of salt restoring to water flux
vprec = forcing.variables['vprec'][n]
# net mass flux crossing surface (including exchange with sea-ice)
PRCmE = forcing.variables['PRCmE'][n]
#--------------------------------------------------------------
# heat flux crossing ocean surface and bottom [Watt/m^2]
# positive values indicate heat entering ocean;
# negative values indicate heat leaving ocean.
# geothermal heating at ocean bottom
#geothermal = forcing.variables['internal_heat'][n]
# net heat crossing ocean surface due to all processes, except restoring
net_heat_surface = forcing.variables['net_heat_surface'][n]
# net heat passed through coupler from shortwave, longwave, latent, sensible.
# note: latent includes heat to vaporize liquid and heat to melt ice/snow.
# note: sensible includes air-sea and ice-sea sensible heat fluxes.
net_heat_coupler = forcing.variables['net_heat_coupler'][n]
# sum of longwave + latent + sensible
LwLatSens = forcing.variables['LwLatSens'][n]
# net shortwave passing through ocean surface
SW = forcing.variables['SW'][n]
# heating of liquid seawater due to formation of frazil sea ice
frazil = forcing.variables['frazil'][n]
# heat from melting/forming sea ice
melth = forcing.variables['seaice_melt_heat'][n]
# there is no restoring heat flux with this test case
heat_restore = 0.0*forcing.variables['frazil'][n]
# net heat content associated with transfer of mass across ocean surface,
# computed relative to 0C. Both diagnostics should be the same, though
# they are computed differently in MOM6.
heat_pme = forcing.variables['Heat_PmE'][n]
heat_content_surfwater = forcing.variables['heat_content_surfwater'][n]
# heat content associated with water mass leaving ocean
heat_content_massout = forcing.variables['heat_content_massout'][n]
# heat content associated with water mass entering ocean
heat_content_massin = forcing.variables['heat_content_massin'][n]
# heat content associated with liquid precipitation
heat_content_lprec = forcing.variables['heat_content_lprec'][n]
# heat content associated with evaporation
heat_content_evap = forcing.variables['heat_content_evap'][n]
# heat content associated with meltw
# no longer needed (GMM)
#heat_content_meltw = forcing.variables['heat_content_icemelt'][n]
# heat content associated with frozen precipitation
heat_content_fprec = forcing.variables['heat_content_fprec'][n]
# heat content associated with virtual precipitation
heat_content_vprec = forcing.variables['heat_content_vprec'][n]
# heat content associated with liquid runoff
heat_content_lrunoff = forcing.variables['heat_content_lrunoff'][n]
# heat content associated with frozen runoff
heat_content_frunoff = forcing.variables['heat_content_frunoff'][n]
# heat content associated with liquid runoff
heat_content_frunoff_glc = forcing.variables['heat_content_frunoff_glc'][n]
# heat content associated with frozen runoff
heat_content_lrunoff_glc = forcing.variables['heat_content_lrunoff_glc'][n]
# heat content associated with liquid condensation
heat_content_cond = forcing.variables['heat_content_cond'][n]
#--------------------------------------------------------------
# salt flux crossing ocean surface and bottom [kg/(m^2 s)]
# positive values indicate salt entering ocean;
# negative values indicate salt leaving ocean.
salt_flux = forcing.variables['salt_flux_in'][n]
# salt flux associated with surface restoring.
# salt_flux has contribution from sea ice + restoring, so we need to remove salt_flux (salt_flux_in)
salt_restore = forcing.variables['salt_flux'][n] - salt_flux# for easy setup of subplots
def newSP(y,x):
global __spv, __spi ; __spv = (y,x) ; __spi = 1 ; plt.subplot(__spv[0], __spv[1], __spi)
def nextSP():
global __spv, __spi ; __spi = __spi + 1 ; plt.subplot(__spv[0], __spv[1], __spi)
# to reduce the amount of code when plotting fields
def make_plot(lon, lat, field, title, xmin=-280, xmax=70, ymin=-80, ymax=90, cmin=-200, cmax=200, xlabel=False):
'''
Uses pcolormesh to plot field as a function of lat and lon.
Writes the max, min and ave values of field into the plot.
'''
global area
field_min = numpy.amin(field)
field_max = numpy.amax(field)
field_ave = (field*area).sum() / area.sum()
ch = plt.pcolormesh(lon,lat,field)
cbax=plt.colorbar(ch, extend='both')
plt.title(r''+title)
if (cmin != 0.0 or cmax != 0.0):
plt.clim(cmin,cmax)
plt.xlim(xmin,xmax)
plt.ylim(ymin,ymax)
plt.ylabel(r'Latitude [$\degree$N]')
if xlabel: plt.xlabel(r'Longitude')
axis = plt.gca()
axis.annotate('max=%5.2f\nmin=%5.2f\nave=%5.2f'%(field_max,field_min,field_ave),xy=(0.01,0.74),
xycoords='axes fraction', verticalalignment='bottom', fontsize=8, color='black')Mass fluxes and global seawater mass budget
Global seawater mass budget consistency check
We compute the change in seawater mass over a given time period. Two different methods are used, and the two methods should agree at the level of truncation error. Note that "truncation error" precision is somewhat larger using offline diagnostics relative to online calculations, particularly if single precision output is saved rather than double precision.
The net mass per time of water (units of kg/s) entering through the ocean boundaries is given by the area integral boundary water mass entering liquid seawater = ∫QW dA, where the net water flux (units of kg m−2 s−1) is given by $$\begin{align*} Q_{W} &= {\tt PRCmE} \end{align*}$$ A nonzero surface mass flux is associated with liquid and solid precipitation and runoff; evaporation and condensation; sea ice melt/formation; and surface restoring.
The time change of liquid seawater mass is computed according to
$$\begin{equation*}
\mbox{seawater mass change} =
\frac{1}{\tau_{n+1} - \tau_{n} } \int dA \left(\int (\rho_{n+1} -
\rho_{n}) \, \mathrm{d}z \right)
\end{equation*}$$ where τn + 1 − τn
is the time increment in seconds. Note that we make use of the MOM6
diagnostic for depth integrated density $$\begin{equation*}
{\tt mass\_wt} = \int \rho \, \mathrm{d}z.
\end{equation*}$$ For a Boussinesq fluid, the in-situ ρ factor is set to ρ0, in which case the
diagnostic field {\tt mass_wt} measures the thickness of a fluid column,
multiplied by ρ0.
For self-consistency, we should have the following equality holding to
within truncation error
$$\begin{equation*}
\boxed{
\mbox{boundary water mass entering liquid seawater} = \mbox{seawater
mass change}.
}
\end{equation*}$$
n0 = n-1
dmass_wt = mass_wt[n] - mass_wt[n0]
dt = time[n] - time[n0]
lhs = area * dmass_wt / dt
rhs = area * ( PRCmE )
print ('Total seawater mass at time step n [kg seawater] =',(mass_wt[n]*area).sum())
print ('Total seawater mass at time step n0 [kg seawater] =',(mass_wt[n0]*area).sum())
print ('Total seawater mass content change [kg seawater] =',dt*lhs.sum())
print ('Net water mass through boundaries [kg seawater] =',dt*rhs.sum())
print ('Residual [kg seawater] =',dt*lhs.sum() - dt*rhs.sum())
print ('Non-dimensional residual (based on difference) =',( lhs.sum() - rhs.sum() )/lhs.sum())
print ('Non-dimensional residual (based on mass_wt[n]) =',( lhs.sum() - rhs.sum() )/(mass_wt[n]*area).sum())Total seawater mass at time step n [kg seawater] = 1.3839304418074236e+21
Total seawater mass at time step n0 [kg seawater] = 1.3839298418176187e+21
Total seawater mass content change [kg seawater] = 599989804701199.2
Net water mass through boundaries [kg seawater] = 599989804700431.5
Residual [kg seawater] = 767.75
Non-dimensional residual (based on difference) = 1.2795169810210496e-12
Non-dimensional residual (based on mass_wt[n]) = 6.4203975413335464e-24
Surface mass fluxes I: combined fields
plt.figure(figsize=(18,8))
newSP(2,2);
field = 86400.0*PRCmE
make_plot(lon,lat,field, '$PRCmE$ [$kg/m^2/day$]',cmin=-20,cmax=20)
nextSP()
field = 86400.0*net_massin
make_plot(lon,lat,field, 'net_massin [$kg/m^2/day$]',cmin=-20,cmax=20)
nextSP()
field = 86400.0*net_massout
make_plot(lon,lat,field, 'net_massout [$kg/m^2/day$]',cmin=-20,cmax=20,xlabel=True)
nextSP()
field = 86400.0*(PRCmE - net_massout - net_massin)
make_plot(lon,lat,field, '$PRCmE - M_{in} - M_{out}$ [$kg/m^2/day$]',cmin=-1e-13,cmax=1e-13,xlabel=True)
Surface mass fluxes II: component fields
plt.figure(figsize=(20,9.5))
newSP(3,3);
field = 86400.0*frunoff_glc
make_plot(lon,lat,field, 'frunoff_glc [$kg/m^2/day$]',cmin=0,cmax=0)
nextSP()
field = 86400.0*lrunoff_glc
make_plot(lon,lat,field, 'lrunoff_glc [$kg/m^2/day$]',cmin=-1.0e-5,cmax=1.0e-5)
nextSP()
field = 86400.0*frunoff
make_plot(lon,lat,field, 'frunoff [$kg/m^2/day$]',cmin=-1,cmax=1)
nextSP()
field = 86400.0*lrunoff
make_plot(lon,lat,field, 'lrunoff [$kg/m^2/day$]',cmin=-10,cmax=10)
nextSP()
field = 86400.0*meltw
make_plot(lon,lat,field, 'Meltw [$kg/m^2/day$]',cmin=-10,cmax=10)
nextSP()
# lprec contains a contribution from sea ice melt/formation.
field = 86400.0*lprec
make_plot(lon,lat,field, 'lprec [$kg/m^2/day$]',cmin=-20,cmax=20)
nextSP()
field = 86400.0*fprec
make_plot(lon,lat,field, 'fprec [$kg/m^2/day$]',cmin=-1,cmax=1,xlabel=True)
nextSP()
# evaporation and condensation
field = 86400.0*evap
make_plot(lon,lat,field, 'evap [$kg/m^2/day$]',cmin=-20,cmax=20,xlabel=True)
nextSP()
field = 86400.0*vprec
make_plot(lon,lat,field, 'vprec [$kg/m^2/day$]',cmin=-5,cmax=5,xlabel=True)
Surface mass flux self-consistency check
We will now check if PRCmE = lprec + fprec + lrunoff + frunoff + vprec + evap + meltw + frunoff_glc + lrunoff_glc
plt.figure(figsize=(10,9))
newSP(2,1);
# this should be within "truncation error" precision
field = 86400.0*(PRCmE -lprec -fprec -lrunoff -frunoff -vprec -evap - meltw - frunoff_glc - lrunoff_glc)
make_plot(lon,lat,field, 'PRCmE -lprec -fprec -lrunoff -frunoff -vprec -evap -meltw [$kg/m^2/day$]',cmin=0.0,cmax=0.0)
Heat fluxes and global ocean heat budget
Global heat budget consistency check
We compute the change in seawater heat content over a given time period. Two different methods are used, and the two methods should agree at the level of truncation error. If larger differences exist, then there is a bug.
The net heat per time (units of Watts) entering through the ocean
boundaries is given by the area integral boundary heating of liquid
seawater = ∫Q dA, where the net heat
flux (units of W m−2 s−1) is given by
$$\begin{align*}
Q &= {\tt (net\_heat\_coupler + heat\_pme + frazil) + internal\_heat
+ heat\_restore} \\
&= {\tt net\_heat\_surface + internal\_heat + heat\_restore}
\end{align*}$$ The time change of liquid seawater heat is
computed according to $$\begin{equation*}
\mbox{seawater heat content change} =
\frac{C_p }{\tau_{n+1} - \tau_{n} } \int dA \left(\rho_0 \int
(\Theta_{n+1} - \Theta_{n}) \, \mathrm{d}z \right)
\end{equation*}$$ where τn + 1 − τn
is the time increment in seconds. Note that we make use of the MOM6
diagnostic for depth integrated potential/conservative temperature
$$\begin{equation*}
{\tt tomint} = \rho_0 \int \Theta \, \mathrm{d}z,
\end{equation*}$$ where the Boussinesq reference density, ρ0, is used since this
test case makes the Boussinesq approximation. For self-consistency, we
should have the following equality holding to within truncation
error
$$\begin{equation*}
\boxed{
\mbox{boundary heating of liquid seawater} = \mbox{seawater heat content
change}.
}
\end{equation*}$$
n0 = n-1
dtomint = tomint[n] - tomint[n0]
dt = time[n] - time[n0]
lhs = cp * area * dtomint / dt
rhs = area * ( net_heat_coupler + heat_pme + frazil)
print ('Total seawater heat at time step n [Joules] =',cp * (area * tomint[n]).sum())
print ('Total seawater heat at time step n0 [Joules] =',cp * (area* tomint[n0]).sum())
print ('Total seawater heat content change [Joules] =',dt*lhs.sum())
print ('Net heat through boundaries [Joules] =',dt*rhs.sum())
print ('Residual [Joules] =',dt*lhs.sum() - dt*rhs.sum())
print ('Non-dimensional residual (based on difference) =',( lhs.sum() - rhs.sum() )/lhs.sum())
print ('Non-dimensional residual (based on tomint[n]) =',( lhs.sum() - rhs.sum()) / (cp * (area * tomint[n]).sum()))Total seawater heat at time step n [Joules] = 1.991087378480628e+25
Total seawater heat at time step n0 [Joules] = 1.9910729648317823e+25
Total seawater heat content change [Joules] = 1.4413648845979715e+20
Net heat through boundaries [Joules] = 1.4118478295111918e+20
Residual [Joules] = 2.951705508677976e+18
Non-dimensional residual (based on difference) = 0.02047854460878781
Non-dimensional residual (based on tomint[n]) = 1.7158090886177494e-12
Basic components to surface heat flux
Self-consistency check: net_heat_surface = heat_pme + frazil + net_heat_coupler
plt.figure(figsize=(18,12))
newSP(3,2);
field = net_heat_surface
make_plot(lon,lat,field, 'net_heat_surface[$W/m^2$]')
nextSP()
field = net_heat_coupler
make_plot(lon,lat,field, 'net_heat_coupler[$W/m^2$]')
nextSP()
field = frazil
make_plot(lon,lat,field, 'frazil[$W/m^2$]',cmin=-20, cmax=20)
nextSP()
field = heat_pme
make_plot(lon,lat,field, 'heat_pme[$W/m^2$]',cmin=-20,cmax=20,xlabel=True)
nextSP()
field = net_heat_surface-net_heat_coupler-frazil-heat_pme
make_plot(lon,lat,field, 'Residual(error)[$W/m^2$]',cmin=0.0,cmax=0.0,xlabel=True)
Heat fluxes crossing ocean surface via the coupler
Heat fluxes crossing ocean surface via the coupler: net_heat_coupler = LwLatSens + SW + melth, where LwLatSens = LW + Latent + Sensible
plt.figure(figsize=(16,12))
newSP(3,2);
field = net_heat_coupler
make_plot(lon,lat,field, 'net_heat_coupler[$W/m^2$]',cmin=-200,cmax=200)
nextSP()
field = LwLatSens
make_plot(lon,lat,field, 'LwLatSens [$W/m^2$]',cmin=-200,cmax=200)
nextSP()
field = SW
make_plot(lon,lat,field, 'SW [$W/m^2$]',cmin=-200,cmax=200)
nextSP()
field = melth
make_plot(lon,lat,field, 'Melth [$W/m^2$]',cmin=-200,cmax=20, xlabel=True)
nextSP()
field = net_heat_coupler - SW - LwLatSens - melth
make_plot(lon,lat,field, 'Residual(error) [$W/m^2$]',cmin=0.0,cmax=0.0, xlabel=True)
Relation between heat_PmE, heat_content_massin, and heat_content_massout
Alternative means to compute to heat_PmE via heat_content_massin and heat_content_massout, where heat_PmE = heat_content_massin + heat_content_massout
This is no longer valid in CESM/MOM6 when ENTHALPY_FROM_COUPLER = True. <div class="alert-warning>
plt.figure(figsize=(18,12))
newSP(3,2);
field = heat_content_massout
make_plot(lon,lat,field, 'heat_content_massout [$W/m^2$]',cmin=-20,cmax=0)
nextSP()
field = heat_content_massin
make_plot(lon,lat,field, 'heat_content_massin [$W/m^2$]',cmin=0,cmax=20)
nextSP()
field = heat_content_massin + heat_content_massout
make_plot(lon,lat,field, 'heat_content_massin + heat_content_massout [$W/m^2$]',cmin=-20,cmax=20)
nextSP()
field = heat_pme
make_plot(lon,lat,field, 'heat_pme [$W/m^2$]',cmin=-20,cmax=20,xlabel=True)
nextSP()
# this should be within "truncation error" precision
field = heat_content_massout + heat_content_massin - heat_pme
make_plot(lon,lat,field, 'heat_massin + heat_massout - heat_pme [$W/m^2$]',cmin=0.0,cmax=0.0, xlabel=True)/glade/derecho/scratch/gmarques/tmp/ipykernel_19633/338529236.py:29: UserWarning: Warning: converting a masked element to nan.
axis.annotate('max=%5.2f\nmin=%5.2f\nave=%5.2f'%(field_max,field_min,field_ave),xy=(0.01,0.74),
Components of heat content from surface mass fluxes
Components of heat content of surface mass fluxes: heat_PmE = heat_content_lprec + heat_content_fprec + heat_content_vprec + heat_content_lrunoff + heat_content_frunoff +heat_content_meltw + heat_content_cond + heat_content_massout + heat_content_frunoff_glc + heat_content_lrunoff_glc
plt.figure(figsize=(20,10))
newSP(3,3);
field = heat_content_frunoff_glc
make_plot(lon,lat,field, 'heat_content_frunoff_glc [$W/m^2$]',cmin=-0.0,cmax=0.0)
nextSP()
field = heat_content_lrunoff_glc
make_plot(lon,lat,field, 'heat_content_lrunoff_glc [$W/m^2$]',cmin=0.0,cmax=0.0)
nextSP()
field = heat_content_lprec
make_plot(lon,lat,field, 'heat_content_lprec [$W/m^2$]',cmin=-20.0,cmax=20.0)
nextSP()
field = heat_content_lrunoff
make_plot(lon,lat,field, 'heat_content_lrunoff [$W/m^2$]',cmin=-20.0,cmax=20.0)
nextSP()
field = heat_content_frunoff
make_plot(lon,lat,field, 'heat_content_frunoff [$W/m^2$]',cmin=0.0,cmax=0.0)
nextSP()
field = heat_content_cond
make_plot(lon,lat,field, 'heat_content_cond [$W/m^2$]',cmin=-1.0,cmax=1.0)
nextSP()
field = heat_content_fprec
make_plot(lon,lat,field, 'heat_content_fprec [$W/m^2$]',cmin=-0.2,cmax=0.2,xlabel=True)
nextSP()
field = heat_content_vprec
make_plot(lon,lat,field, 'heat_content_vprec [$W/m^2$]',cmin=-0.0,cmax=.0,xlabel=True)
Self-consistency of diagnosed heat content from mass entering ocean
This is no longer valid in CESM/MOM6 when ENTHALPY_FROM_COUPLER = True. <div class="alert-warning>
plt.figure(figsize=(16,12))
newSP(3,2);
heat_content_sum = ( heat_content_lprec + heat_content_fprec + heat_content_vprec +
heat_content_lrunoff + heat_content_cond + heat_content_frunoff +
heat_content_frunoff_glc + heat_content_lrunoff_glc)
field = heat_content_massin
make_plot(lon,lat,field, 'heat_content_massin [$W/m^2$]',cmin=-20.0,cmax=20.0)
nextSP()
field = heat_content_sum
make_plot(lon,lat,field, '$\Sigma$ of heat_contents entering ocean [$W/m^2$]',cmin=-20.0,cmax=20.0,xlabel=True)
nextSP()
# this should be within "truncation error" precision
field = heat_content_massin - heat_content_sum
make_plot(lon,lat,field, 'heat_content_massin - $\Sigma$ of components [$W/m^2$]',cmin=0.0,cmax=0.0,xlabel=True)/glade/derecho/scratch/gmarques/tmp/ipykernel_19633/338529236.py:29: UserWarning: Warning: converting a masked element to nan.
axis.annotate('max=%5.2f\nmin=%5.2f\nave=%5.2f'%(field_max,field_min,field_ave),xy=(0.01,0.74),
Self-consistency between heat_pme and heat_content_surfwater
comp_sum = ( heat_content_lprec + heat_content_fprec + heat_content_vprec +
+ heat_content_frunoff + heat_content_cond + heat_content_evap +
heat_content_lrunoff + heat_content_lrunoff_glc + heat_content_frunoff_glc)
plt.figure(figsize=(16,12))
newSP(3,2);
field = heat_content_surfwater
make_plot(lon,lat,field, 'heat_content_surfwater [$W/m^2$]',cmin=-20.0,cmax=20.0)
nextSP()
field = comp_sum
make_plot(lon,lat,field, '$\Sigma$ of heat_content components [$W/m^2$]',cmin=-20.0,cmax=20.0)
nextSP()
# this should be within "truncation error" precision
field = comp_sum - heat_content_surfwater
make_plot(lon,lat,field, '$\Sigma$ - heat_content_surfwater [$W/m^2$]',cmin=0.0,cmax=0.0)
nextSP()
field = heat_pme
make_plot(lon,lat,field, 'heat_pme [$W/m^2$]',cmin=-20.0,cmax=20.0, xlabel=True)
nextSP()
# this should be within "truncation error" precision
field = heat_pme - heat_content_surfwater
make_plot(lon,lat,field, 'heat_pme - heat_content_surfwater [$W/m^2$]',cmin=0.0,cmax=0.0, xlabel=True)
Map effective temperatures
The following "effective" temperatures differ generally from the SST due to the means by which the water is exchanged across the ocean surface boundary. In particular, there are some occasions whereby layers deeper than k=1 need to be sampled in order to exchange water with the atmosphere and/or sea ice. These maps provide us with a sense for where such occurrances take place.
This is no longer valid in CESM/MOM6 when ENTHALPY_FROM_COUPLER = True. <div class="alert-warning>
TinEff = heat_content_massin/(net_massin*cp)
ToutEff = heat_content_massout/(net_massout*cp)
TnetEff = heat_pme/(PRCmE*cp)
plt.figure(figsize=(16,12))
newSP(3,2);
field = TinEff
make_plot(lon,lat,field, '$\Theta_{in} [{\degree}C$]',cmin=-2,cmax=30)
nextSP()
field = TinEff - sst
make_plot(lon,lat,field, '$\Theta_{in} - SST [{\degree}C$]',cmin=0.0,cmax=0.0)
nextSP()
field = ToutEff
make_plot(lon,lat,field, '$\Theta_{out} [{\degree}C$]',cmin=-2,cmax=30)
nextSP()
field = ToutEff - sst
make_plot(lon,lat,field, '$\Theta_{out} - SST [{\degree}C$]',cmin=0.0,cmax=0.0)
nextSP()
field = TnetEff
make_plot(lon,lat,field, '$\Theta_{net} [{\degree}C$]',cmin=-2,cmax=30, xlabel=True)
nextSP()
field = TnetEff - sst
make_plot(lon,lat,field, '$\Theta_{net} - SST [{\degree}C$]',cmin=0.0,cmax=0.0, xlabel=True)/glade/derecho/scratch/gmarques/tmp/ipykernel_19633/338529236.py:29: UserWarning: Warning: converting a masked element to nan.
axis.annotate('max=%5.2f\nmin=%5.2f\nave=%5.2f'%(field_max,field_min,field_ave),xy=(0.01,0.74),
Salt fluxes and global ocean salt budget
Global salt budget consistency check
We compute the change in seawater salt content over a given time period. Two different methods are used, and the two methods should agree at the level of truncation error. If larger differences exist, then there is a bug.
The net salt per time (units of kg/s) entering through the ocean boundaries is given by the area integral boundary salt entering liquid seawater = ∫QS dA, where the net salt flux (units of kg m−2 s−1) is given by $$\begin{align*} Q_{S} &= {\tt salt\_flux} + {\tt salt\_restore}. \end{align*}$$ A nonzero salt flux is associated with exchanges between liquid seawater and solid sea ice. It also arises from simulations using a restoring boundary flux associated with damping to observed sea surface salinity. Finally, there can be a salt flux when using sponges to damp the ocean interior back to an observed value.
The time change of liquid seawater salt content is computed according
to $$\begin{equation*}
\mbox{seawater salt content change} =
\frac{1}{\tau_{n+1} - \tau_{n} } \int dA \left(\rho_0 \int (S_{n+1} -
S_{n}) \, \mathrm{d}z \right)
\end{equation*}$$ where τn + 1 − τn
is the time increment in seconds. Note that we make use of the MOM6
diagnostic for depth integrated salinity $$\begin{equation*}
{\tt somint} = \rho_0 \int S \, \mathrm{d}z,
\end{equation*}$$ where the Boussinesq reference density, ρ0, is used since this
test case makes the Boussinesq approximation. For self-consistency, we
should have the following equality holding to within truncation
error
$$\begin{equation*}
\boxed{
\mbox{boundary salt entering liquid seawater} = \mbox{seawater salt
content change}.
}
\end{equation*}$$
n0 = n-1
dsomint = somint[n] - somint[n0]
time = surface.variables[tvar][:]*86400.
dt = time[n] - time[n0]
lhs = 1.e-3 * area * dsomint / dt
rhs = area * ( salt_flux + salt_restore )
print ('Total seawater salt at time step n [kg salt] =',(area * somint[n]).sum())
print ('Total seawater salt at time step n0 [kg salt] =',(area * somint[n0]).sum())
print ('Total seawater salt content change [kg salt] =',dt*lhs.sum())
print ('Net salt through boundaries [kg salt] =',dt*rhs.sum())
print ('Residual [kg salt] =',dt*lhs.sum() - dt*rhs.sum())
print ('Non-dimensional residual (based on diference) =',( lhs.sum() - rhs.sum() )/lhs.sum())
print ('Non-dimensional residual (based on somint[n]) =',( lhs.sum() - rhs.sum() )/((area * somint[n]).sum()))Total seawater salt at time step n [kg salt] = 4.805039270307199e+22
Total seawater salt at time step n0 [kg salt] = 4.805039106618365e+22
Total seawater salt content change [kg salt] = 1636888334022.601
Net salt through boundaries [kg salt] = 1636888333948.7808
Residual [kg salt] = 73.8203125
Non-dimensional residual (based on diference) = 4.509800620096413e-11
Non-dimensional residual (based on somint[n]) = 1.778139017120276e-26
plt.figure(figsize=(16,9))
newSP(2,2)
field = 86400.0*salt_flux
make_plot(lon,lat,field, 'Surface salt flux from ice-ocean exchange [kg m$^{-2}$ day$^{-1}$]',cmin=-0.0,cmax=0.0,xlabel=True)
nextSP()
field = 86400.0*(salt_restore)
make_plot(lon,lat,field, 'Surface salt flux from restoring [kg m$^{-2}$ day$^{-1}$]',cmin=-0.0,cmax=0.0,xlabel=True)