FFT-K coefficient question

From: Ventrice, Michael J <mventrice_at_nyahnyahspammersnyahnyah>
Date: Fri Feb 24 2012 - 08:38:08 MST

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