[ncl-talk] Sverdrup stream function
Chathurika
chatu at scsio.ac.cn
Wed Oct 27 08:12:13 MDT 2021
Dear all,
Thank you so much for reading the email.
I am trying to calculate the Sverdrup stream function (using the following equation)
However, My results are not correct. Could you please check the attached script and let me know where I do wrong. That would be a really great help. Please be kind to check this out.
The script:
diri="./"
sfile = "HR_hist_r1i1p1f1.nc"
f = addfile(diri+sfile,"r")
tauuo = f->tauuo(1740:1979,:,:)
printVarSummary(tauuo)
printMinMax(tauuo, False)
tauvo = f->tauvo(1740:1979,:,:)
printVarSummary(tauvo)
printMinMax(tauvo, False)
;***************wind stress curl*********************************
; Miscellaneous
;************************************************
lat = tauuo&lat
lon = tauuo&lon
;************************************************
; CENTERED FINITE DIFFERENCES:
; Use local function [grad_latlon_cfd]
; to compute zonal (X) and meridional (Y) gradients.
;************************************************
gradLatLon = grad_latlon_cfd(tauuo, lat, lon, True, False)
tauuoY_cfd = gradLatLon[0]
tauuoX_cfd = gradLatLon[1]
tauuoY_cfd at long_name = "tauuoY: cfd"
tauuoY_cfd at units = "Nm-3"
tauuoX_cfd at long_name = "tauuoX: cfd"
tauuoX_cfd at units = "Nm-3"
printVarSummary(tauuoY_cfd)
printMinMax(tauuoY_cfd, False) ; unscaled = True
printVarSummary(tauuoX_cfd)
printMinMax(tauuoX_cfd,False)
print("")
delete([/lat,lon,gradLatLon/])
;************************************************
; Miscellaneous
;************************************************
lat = tauvo&lat
lon = tauvo&lon
;************************************************
; CENTERED FINITE DIFFERENCES:
; Use local function [grad_latlon_cfd]
; to compute zonal (X) and meridional (Y) gradients.
;************************************************
gradLatLon = grad_latlon_cfd(tauvo, lat, lon, True, False)
tauvoY_cfd = gradLatLon[0]
tauvoX_cfd = gradLatLon[1]
tauvoY_cfd at long_name = "tauvoY: cfd"
tauvoY_cfd at units = "Nm-3"
tauvoX_cfd at long_name = "tauvoX: cfd"
tauvoX_cfd at units = "Nm-3"
printVarSummary(tauuoY_cfd)
printMinMax(tauvoY_cfd, False) ; True=unscaled
printVarSummary(tauuoX_cfd)
printMinMax(tauvoX_cfd,False)
print("")
curl = (tauvoX_cfd-tauuoY_cfd)
copy_VarCoords(tauuo, curl)
curl at long_name = "wind stress curl"
curl at units = "Nm-3"
printVarSummary(curl)
printMinMax(curl, False)
curl_av = dim_avg_n_Wrap(curl, 0)
printVarSummary(curl_av)
printMinMax(curl_av, False)
delete([/gradLatLon,curl/])
;..........Sv stream function........
re = 6371000
density = 1027
omega = 7.292e-05
dx = 104492.5 ; in meter
d2rad = 0.017453 ; degrees to radians
lat_S = sin(lat*d2rad)
printMinMax(lat_S, True)
coriolis = 2*omega*lat_S
printVarSummary(coriolis)
printMinMax(coriolis, True)
nlat = dimsizes(curl_av&lat)
print(nlat)
nlon = dimsizes(curl_av&lon)
print(nlon)
latitude = curl_av&lat
longitude = curl_av&lon
;**************************************************************
; calculate intergral
; int[lon1:lon2]curl*re*cos(lat)*dx
;**************************************************************
zone_int = new((/nlat,nlon/),typeof(curl_av)) ; allocate space
do k = 0, nlon-1
do j = 0, nlat-1
zone_int(j,k) = dim_sum(curl_av(j,0:k)*re*cos(lat(j)*d2rad)*dx)
end do
end do
printVarSummary(zone_int)
printMinMax(zone_int, True)
cornew = new(dimsizes(curl_av), typeof(coriolis), getFillValue(coriolis))
do i = 0,180
do j = 0,360
cornew(i,j) = coriolis(i)
end do
end do
printVarSummary(cornew)
printMinMax(cornew, True)
dimcornew = dimsizes(cornew)
;************************************************
; SPHERICAL HARMONICS:
; Use "highly accurate" spherical harmonics procedure (gradsg)
; to compute zonal (X) and meridional (Y) gradients.
;************************************************
; pre-allocate space for return gradients
cornewX_gradsg = new( dimcornew, typeof(cornew), getFillValue(cornew) ) ; lon=>X
cornewY_gradsg = new( dimcornew, typeof(cornew), getFillValue(cornew) ) ; lat=>Y
gradsg(cornew, cornewX_gradsg, cornewY_gradsg) ; procedure for gaussian grids
copy_VarCoords(cornew, cornewX_gradsg) ; add meta data
copy_VarCoords(cornew, cornewY_gradsg)
cornewX_gradsg at long_name = "cornewX: gradsg"
cornewX_gradsg at units = "1/s2"
cornewY_gradsg at long_name = "cornewY: gradsg"
cornewY_gradsg at units = "1/s2"
print("")
printMinMax(cornewY_gradsg,True ) ; unscaled
printMinMax(cornewX_gradsg,False)
print("")
printVarSummary(cornewY_gradsg)
;print(cornewY_gradsg)
;
betarow = (cornewY_gradsg*density)
copy_VarCoords(curl_av, betarow)
printVarSummary(betarow)
printMinMax(betarow, True)
stm = (-1/betarow) * zone_int *1e-06
copy_VarCoords(curl_av, stm)
printVarSummary(stm)
printMinMax(stm, True)
;print(stm(80,:))
;**********************plot**************************
wks = gsn_open_wks("png","HR_sv") ; send graphics to PNG file
setvalues NhlGetWorkspaceObjectId
"wsMaximumSize" : 300000000
end setvalues
res = True
res at gsnDraw = False ;-- don't draw plot yet
res at gsnFrame = False ;-- don't advance frame
res at gsnMaximize = True ;-- maximize plot output
res at gsnAddCyclic = False ;-- don't add cyclic points
res at gsnLeftStringFontHeightF = 0.010 ;-- decrease left string font size
res at gsnRightStringFontHeightF = 0.010 ;-- decrease right string font size
res at cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
res at cnMinLevelValF = -80 ; set min contour level
res at cnMaxLevelValF = 80 ; set max contour level
res at cnLevelSpacingF = 1 ; set contour spacing
res at cnFillOn = True ; turn on contour fill
res at cnLinesOn = False ; turn off contour lines
res at cnLineLabelsOn = False ; turn off contour labels
res at cnFillDrawOrder = "PreDraw"
res at cnFillPalette = "cmp_flux" ; set color map
cmap = read_colormap_file("cmp_flux") ; 16 x 4
ncolors = dimsizes(cmap(:,1))
print(dimsizes(cmap))
cmap(10,:) = namedcolor2rgba("white")
cmap(11,:) = namedcolor2rgba("white")
res at cnFillPalette := cmap
res at mpFillOn = True
;res at mpOutlineOn=False
res at mpGridAndLimbOn = False
res at mpOutlineBoundarySets = "National"
res at mpGeophysicalLineThicknessF = 1.5
res at mpDataBaseVersion = "LowRes"
res at mpDataSetName = "Earth.4"
res at lbLabelBarOn = True
;res at sfXArray = lon2d ;-- longitude grid cell center
;res at sfYArray = lat2d ;-- latitude grid cell center
;-- create the plot
res at gsnCenterString = "Wind stress curl"
res at gsnCenterStringFontHeightF = 0.013
res at gsnLeftStringFontHeightF = 0.013
res at gsnLeftString = "[2080-2099]-[1995-2014]"
res at gsnRightStringFontHeightF = 0.013
res at gsnRightString = "MPI-ESM-HR"
plot = gsn_csm_contour_map(wks, stm, res)
;-- draw the plot and advance the frame
draw(plot)
frame(wks)
Thank you very much and best regards,
Chathu
Wickramage Chathurika Hemamali
Msc in Physical Oceanography
State Key Laboratory of Tropical Oceanography
South China Sea Institute of Oceanology
University of Chinese Academy of Science
China
Email : wickramagechathurika at rocketmail.com
chatu at scsio.ac.cn
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20211027/cd8cd540/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 1635343359227.png
Type: image/png
Size: 3034 bytes
Desc: not available
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20211027/cd8cd540/attachment.png>
More information about the ncl-talk
mailing list