Hello,
I am trying to initialize a barotropic model using a set of initial conditions. I can initialize the vorticity no problem, however when I try to calculate the stream function, I am having problems understanding how to set up my K-coefficient array. Could you please help me try to figure out how to set up my K-array?
Thank you very much,
-Mike
My Script:
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
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"
begin
plotType = "x11"
plotName = "sf"
fontHeightF = 0.02
; ##### Set initial conditions for barotropic model #####
Lx = 10.
Ly = 10.
Nx = 64
Ny = 64
Nsteps = 64
Tmax = 60
PI = 3.14159
; ##### to make X and Y 1d arrays #####
X1 = fspan(-Lx/2, Lx/2, Nx)
Y1 = fspan(-Ly/2, Ly/2, Ny)
X_ar = new((/Nx,Ny/), "float")
X = conform_dims(dimsizes(X_ar), X1, 1)
Y_ar = new((/Nx,Ny/), "float")
Y = conform_dims(dimsizes(Y_ar), Y1, 0 )
; ##### Create Vorticity #####
Zeta = new((/Tmax, Nx, Ny/), "float")
psi = new((/Tmax, Nx, Ny/), "float")
Zeta(0,:,:) = exp(-((X^2)/4) - Y^2) ;Initial time
; ##### Calculate wave number coefficients #####
; ##### where k = 2PI/L #####
L = 10
K = new((/Nx,Ny/2+1/), "float")
piLx = 2 * PI / Lx
piLy = 2 * PI / Ly
noZero = 1.0e-7 ;If you change this to 1.0e^7, you get the correct shape for SF, but wrong amplitude
do i = 0, Nx-1
do j = 0, Ny/2
K(i,j) = 1.0 / max((/ (piLx * int2flt(i-Nx/2))^2 \\
+ (piLy * int2flt(j))^2, noZero /))
end do
end do
dy = Ly/Ny
dx = Lx/Nx
dt = 0.375
t = 0 ; Just for initial time
fftcoef = fft2df(Zeta(0,:,:))
do i = 0, Nx-1
do j = 0, Ny/2
fftcoef(:,i,j) = -fftcoef(:,i,j) * K(i,j)
end do
end do
psi(t,:,:) = fft2db(fftcoef)
;--------Resources------------
if( ( plotType.eq."png" ).or.( plotType.eq."gif" ) ) then
plotTypeLocal = "eps"
else
plotTypeLocal = plotType
end if
wks = gsn_open_wks( plotTypeLocal, plotName )
gsn_define_colormap(wks,"amwg_blueyellowred")
res = True
res@cnFillOn<mailto:res@cnFillOn> = True
res@cnMonoFillColor<mailto:res@cnMonoFillColor> = False
res@cnLineLabelsOn<mailto:res@cnLineLabelsOn> = False
res@cnInfoLabelOn<mailto:res@cnInfoLabelOn> = False
res@cnLinesOn<mailto:res@cnLinesOn> = False
res@mpFillOn<mailto:res@mpFillOn> = False
res@cnFillDrawOrder<mailto:res@cnFillDrawOrder> = "PreDraw"
res@mpDataBaseVersion<mailto:res@mpDataBaseVersion> = "MediumRes"
res@tmXBLabelFontHeightF<mailto:res@tmXBLabelFontHeightF> = fontHeightF
res@tmYLLabelFontHeightF<mailto:res@tmYLLabelFontHeightF> = fontHeightF
res@gsnLeftStringFontHeightF<mailto:res@gsnLeftStringFontHeightF> = fontHeightF
res@gsnRightStringFontHeightF<mailto:res@gsnRightStringFontHeightF> = fontHeightF
res@lbTitleFontHeightF<mailto:res@lbTitleFontHeightF> = fontHeightF
res@lbLabelFontHeightF<mailto:res@lbLabelFontHeightF> = fontHeightF
res@cnInfoLabelFontHeightF<mailto:res@cnInfoLabelFontHeightF> = fontHeightF
res@lbBoxLinesOn<mailto:res@lbBoxLinesOn> = True
res@lbLabelBarOn<mailto:res@lbLabelBarOn> = True
res@gsnLeftString<mailto:res@gsnLeftString> = "Barotropic Model"
res@gsnRightString<mailto:res@gsnRightString> = "Mike Ventrice"
res@gsnSpreadColors<mailto:res@gsnSpreadColors> = True
res@cnLevelSelectionMode<mailto:res@cnLevelSelectionMode> = "ManualLevels"
res@cnMinLevelValF<mailto:res@cnMinLevelValF> = 0.1
res@cnMaxLevelValF<mailto:res@cnMaxLevelValF> = 1.
res@cnLevelSpacingF<mailto:res@cnLevelSpacingF> = 0.05
resPSI = res
resPSI@cnFillOn<mailto:resPSI@cnFillOn> = True
resPSI@cnMonoFillColor<mailto:resPSI@cnMonoFillColor> = False
resPSI@cnLineLabelsOn<mailto:resPSI@cnLineLabelsOn> = True
resPSI@cnInfoLabelOn<mailto:resPSI@cnInfoLabelOn> = False
resPSI@gsnLeftString<mailto:resPSI@gsnLeftString> = " "
resPSI@gsnRightString<mailto:resPSI@gsnRightString> = " "
resPSI@cnLinesOn<mailto:resPSI@cnLinesOn> = True
resPSI@cnLevelSelectionMode<mailto:resPSI@cnLevelSelectionMode> = "AutomaticLevels"
resPSI@cnLevels<mailto:resPSI@cnLevels> = (/-5,-0.5, -0.05, -0.005, -0.0005, -0.00005/)
resPSI@cnMinLevelValF<mailto:resPSI@cnMinLevelValF> = -0.5
resPSI@cnMaxLevelValF<mailto:resPSI@cnMaxLevelValF> = -0.05
resPSI@cnLevelSpacingF<mailto:resPSI@cnLevelSpacingF> = 0.05
resPSI@gsnSpreadColors<mailto:resPSI@gsnSpreadColors> = True
plotPsi = gsn_csm_contour( wks, psi(0,:,:), resPSI )
system( "date" )
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Feb 24 08:38:29 2012
This archive was generated by hypermail 2.1.8 : Tue Mar 20 2012 - 15:27:15 MDT