# [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"

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 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("")

;************************************************
; 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 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)

;..........Sv stream function........

re = 6371000
density = 1027
omega = 7.292e-05

dx = 104492.5                 ; in meter
d2rad = 0.017453         ; degrees to radians

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
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

cornewX_gradsg at units     = "1/s2"
cornewY_gradsg at units     = "1/s2"

print("")
printMinMax(cornewY_gradsg,True )          ; unscaled
print("")

;

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>
```