NCL Home >
Documentation >
Functions >
Oceanography
rho_mwjf
Computes ocean water density given a specified range for potential temperature (deg Celsius) and salinity (psu).
Prototype
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" ; This library is automatically loaded
; from NCL V6.2.0 onward.
; No need for user to explicitly load.
function rho_mwjf (
t2d [*][*] : numeric,
s2d [*][*] : numeric,
depth [1] : numeric
)
return_val : typeof(t2d)
Arguments
t2dA two-dimensional array of temperature values (must be deg C).
s2dA two-dimensional array of the same size as t2d containing salinity values (must be psu).
depthA scalar value for depth at which to compute density (meters).
Return value
The results are returned in an array of the same type and size as t2d.
Description
This function computes potential density for a specifed range or potential T and Salt. It can also be used for water mass T-S diagrams (see the example below).
This function is based on Steve Yeager's "Rhoalphabeta", which uses the CESM3 POP equation of state "mwjf" calculation to compute density.
Note: In September 2009, Dr. Arne Melson [Norwegian Meteorological Institute, R & D - Oceanography ] noted that the results were applicable for computing surface potential density only, i.e., where depth=0. A modification has since been applied [version 5.2.0] to accurately compute all density surfaces (depth >= 0).
References:
Levitus, S., R. Burgett, and T.P. Boyer, World Ocean Atlas
1994, Volume 3: Salinity, NOAA Atlas NESDIS 3, US Dept. of
Commerce, 1994.
Levitus, S. and T.P. Boyer, World Ocean Atlas 1994,
Volume 4: Temperature, NOAA Atlas NESDIS 4, US Dept. of
Commerce, 1994.
Dukowicz, J. K., 2000: Reduction of Pressure and Pressure
Gradient Errors in Ocean Simulations, J. Phys. Oceanogr.,
submitted.
Examples
; read potential temp (TEMP), salinity (SALT)
; compute potential density (PD) for specified range PD(t,s)
; (use ncl function based yeager's algorithm for rho computation)
; assumes annual and zonally avgeraged input data set (i.e, one time slice)
; used k.lindsay's "za" for zonal avg -- already binned into basins
; plots temp vs salt (scatter plot), pd overlay
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
; ================================> ; PARAMETERS
case = "PHC2_gx1v3"
ocnfile = "za_PHC2_T_S_gx1v3.nc"
depth_min = 14895.82 ; in cm, depth of first layer to be included
depth_max = 537499.9
;
; plot limits
;
smincn = 32.5
smaxcn = 37.0
tmincn = -2.
tmaxcn = 22.
;
; Choose basin index
;
; 0 = global 1 = southern ocean 2 = pacific 3 = indian 6 = atlantic
; 8 = labrador 9 = GIN 10 = arctic
;
bi = 2
;=====> basin check
if(bi.lt.0.or.bi.gt.10) then
print("basin index "+ bi + " not supported")
exit
end if
if(bi.eq.0) then
basin = "Global"
blab = "global"
end if
if(bi.eq.1) then
basin = "Southern Ocean"
blab = "so"
end if
if(bi.eq.2) then
basin = "Pacific Ocean"
blab = "pacific"
end if
if(bi.eq.3) then
basin = "Indian Ocean"
blab = "indian"
end if
if(bi.eq.6) then
basin = "Atlantic Ocean"
blab = "atlantic"
end if
if(bi.eq.8) then
basin = "Labrador Sea"
blab = "lab"
end if
if(bi.eq.9) then
basin = "GIN Sea"
blab = "gin"
end if
if(bi.eq.10) then
basin = "Arctic Ocean"
blab = "arctic"
end if
;=====> initial resource settings
wks = gsn_open_wks("ps","TS_diagram") ; Open a Postscript file
;===== data
focn = addfile(ocnfile, "r")
salt = focn->SALT(0,:,{depth_min:depth_max},:) ;(basins, z_t, lat_t)
temp = focn->TEMP(0,:,{depth_min:depth_max},:)
;====section out choice basin
temp_ba = temp(bi,:,:)
salt_ba = salt(bi,:,:)
;===== put into scatter array format
tdata_ba = ndtooned(temp_ba)
sdata_ba = ndtooned(salt_ba)
ydata = tdata_ba
xdata = sdata_ba
;============== compute potenial density (PD), using rho_mwjf
;
; for potential density, depth = 0. (i.e. density as if brought to surface)
;
;===========================================================================
; WARNING: T-S diagrams use POTENTIAL DENSITY... if set depth to something
; other then 0, then you will be plotting density contours computed for the
; specified depth layer.
;===========================================================================
depth = 0. ;in meters
tspan = fspan(tmincn,tmaxcn,51)
sspan = fspan(smincn,smaxcn,51)
; the more points the better... using yeager's numbers
t_range = new((/51,51/),float)
s_range = new((/51,51/),float)
t_range = conform(t_range,tspan,0)
s_range = conform(s_range,sspan,1)
pd = rho_mwjf(t_range,s_range,depth)
pd!0 = "temp"
pd!1 = "salt"
pd&temp = tspan
pd&salt = sspan
pd = 1000.*(pd-1.) ; Put into kg/m3 pot den units
printVarSummary(pd)
printVarInfo(pd,"rho_mwjf")
;=================Graphics
;--- scatter plot
res = True
res@gsnDraw = False
res@gsnFrame = False
res@xyMarkLineModes = "Markers"
res@xyMarkers = 16
res@xyMarkerColors = "black"
res@pmLegendDisplayMode = "Never"
res@txFontHeightF = 0.01
res@tiXAxisString = salt@units
res@tiXAxisFontHeightF = 0.02
res@tiYAxisString = temp@units
res@tiYAxisFontHeightF = 0.02
res@trXMinF = smincn
res@trXMaxF = smaxcn
res@trYMinF = tmincn
res@trYMaxF = tmaxcn
res@gsnRightString = depth_min/100. + "-"+depth_max/100. +"m"
res@gsnLeftString = basin
plot = gsn_csm_xy(wks,xdata,ydata,res)
;----- pd overlay
resov = True
resov@gsnDraw = False
resov@gsnFrame = False
resov@cnLevelSelectionMode = "AutomaticLevels"
resov@cnInfoLabelOn = "False"
resov@cnLineLabelPlacementMode = "Constant"
resov@cnLineLabelFontHeightF = ".02"
plotpd = gsn_csm_contour(wks,pd,resov)
overlay(plot,plotpd)
; panel resources
pan = True
pan@gsnMaximize = True
pan@gsnPaperOrientation = "portrait"
pan@txFontHeightF = 0.025
pan@txString = case + " ANN AVG: T-S Diagram "
gsn_panel(wks,plot,(/1,1/),pan)
end