Case: b.e30_a06c.cmeps547¶

caseroot: /glade/work/gmarques/cesm.cases/sfc_fluxes_check/b.e30_a06c.cmeps547

Purposes¶

Evaluate surface flux consistency after incorporating separate runoff fluxes for GLC (see https://github.com/ESCOMP/CMEPS/pull/547).

MOM6 diagnostic boundary fluxes of scalars and their global budgets

This notebook is meant to work with MOM6 run via the CESM framework

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:

  1. Maps of surface boundary fluxes of water, heat, and salt crossing into the liquid seawater in MOM6;
  2. 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:

  1. To provide a starting point to document boundary fluxes of scalar fields;
  2. To teach MOM6 users about the boundary fluxes, their patterns, units, and sign conventions;
  3. To perform self-consistency checks to ensure the model is conserving scalar fields;
  4. 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.

In [1]:
# %pylab inline
import matplotlib.pyplot as plt
import netCDF4
import numpy
In [2]:
path    = '/glade/derecho/scratch/gmarques/b.e30_a06c.cmeps547/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
In [3]:
# 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
In [4]:
# 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 $$\begin{equation*} \mbox{boundary water mass entering liquid seawater} = \int Q_{W} \, dA, \end{equation*}$$ where the net water flux (units of $\mbox{kg}~\mbox{m}^{-2}~\mbox{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 $\tau_{n+1} - \tau_{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 $\rho$ factor is set to $\rho_{0}$, in which case the diagnostic field {\tt mass_wt} measures the thickness of a fluid column, multiplied by $\rho_{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*}$$

In [5]:
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.3839304300782226e+21
Total seawater mass at time step n0 [kg seawater]  = 1.3839298270982617e+21
Total seawater mass content change [kg seawater]   = 602979960582475.8
Net water mass through boundaries [kg seawater]    = 602979960583855.9
Residual [kg seawater]                             = -1380.125
Non-dimensional residual (based on difference)     = -2.2888944258941864e-12
Non-dimensional residual (based on mass_wt[n])     = -1.1542520095393669e-23

Surface mass fluxes I: combined fields

In [6]:
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)
No description has been provided for this image

Surface mass fluxes II: component fields

In [7]:
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=0,cmax=0)

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)
No description has been provided for this image

Surface mass flux self-consistency check

We will now check if PRCmE = lprec + fprec + lrunoff + frunoff + vprec + evap + meltw + frunoff_glc + lrunoff_glc

In [8]:
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)
No description has been provided for this image

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 $$\begin{equation*} \mbox{boundary heating of liquid seawater} = \int Q \, dA, \end{equation*}$$ where the net heat flux (units of $\mbox{W}~\mbox{m}^{-2}~\mbox{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 $\tau_{n+1} - \tau_{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, $\rho_{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*}$$

In [9]:
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.9910882721158677e+25
Total seawater heat at time step n0 [Joules]   = 1.991072947849002e+25
Total seawater heat content change [Joules]    = 1.532426686663185e+20
Net heat through boundaries [Joules]           = 1.5029096315770307e+20
Residual [Joules]                              = 2.951705508615422e+18
Non-dimensional residual (based on difference) = 0.019261642558853268
Non-dimensional residual (based on tomint[n])  = 1.715808318496247e-12

Basic components to surface heat flux

Self-consistency check: net_heat_surface = heat_pme + frazil + net_heat_coupler

In [10]:
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)
No description has been provided for this image

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

In [11]:
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)
No description has been provided for this image

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.
In [12]:
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_42725/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),
No description has been provided for this image

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

In [13]:
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)
No description has been provided for this image

Self-consistency of diagnosed heat content from mass entering ocean

This is no longer valid in CESM/MOM6 when ENTHALPY_FROM_COUPLER = True.
In [14]:
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_42725/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),
No description has been provided for this image

Self-consistency between heat_pme and heat_content_surfwater

In [15]:
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)
No description has been provided for this image

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.
In [16]:
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_42725/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),
No description has been provided for this image

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 $$\begin{equation*} \mbox{boundary salt entering liquid seawater} = \int Q_{S} \, dA, \end{equation*}$$ where the net salt flux (units of $\mbox{kg}~\mbox{m}^{-2}~\mbox{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 $\tau_{n+1} - \tau_{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, $\rho_{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*}$$

In [17]:
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.805039269537945e+22
Total seawater salt at time step n0 [kg salt]    = 4.805039106402767e+22
Total seawater salt content change [kg salt]     = 1631351770378.3267
Net salt through boundaries [kg salt]            = 1631351770120.633
Residual [kg salt]                               = 257.693603515625
Non-dimensional residual (based on diference)    = 1.579632042886278e-10
Non-dimensional residual (based on somint[n])    = 6.207158991196285e-26
In [20]:
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)
No description has been provided for this image
In [ ]: