[ncl-talk] Sverdrup stream function
Dennis Shea
shea at ucar.edu
Fri Oct 29 12:17:54 MDT 2021
I have looked at the script through "calculate integral "
I don't see anything obviously wrong. Perhaps plotting the tauvoX_cfd
and tauvoY_cfd
components might give some insight
curl = (tauvoX_cfd-tauuoY_cfd)
On Wed, Oct 27, 2021 at 8:12 AM Chathurika via ncl-talk <
ncl-talk at mailman.ucar.edu> wrote:
>
> 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
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at mailman.ucar.edu
> List instructions, subscriber options, unsubscribe:
> https://mailman.ucar.edu/mailman/listinfo/ncl-talk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20211029/d611790a/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/20211029/d611790a/attachment.png>
More information about the ncl-talk
mailing list