; Mike Ventrice ; ATM 562 ; Numerical Weather Prediction ; Project 1 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 = "png" 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) ;print(X) 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 the zonal and meridional wave numbers ##### ; ##### where k = 2PI/L ##### L = 10 mult_x_part1 = fspan(0, Nx/2-1, Nx/2) mult_x_part2 = fspan(-Nx/2, -1, Nx/2) mult_x = array_append_record(mult_x_part1, mult_x_part2, 0) mult_y_part1 = fspan(0, Ny/2-1, Ny/2) mult_y_part2 = fspan(-Ny/2, -1, Ny/2) mult_y = array_append_record(mult_y_part1, mult_y_part2, 0) kx = ((2*PI)/L) * mult_x ky = ((2*PI)/L) * mult_y KX_ar = new((/Nx, Ny/), "float") KX = conform_dims(dimsizes(KX_ar), kx, 1) ; print(KX) KX(0,0) = 1 * 10^-6 KY_ar = new((/Nx, Ny/), "float") KY = conform_dims(dimsizes(KY_ar), ky, 0) KY(0,0) = 1 * 10^-6 print(KY) K2 = 1/((KX^2)+(KY^2)) ;print(K2) ; print("Ly = "+Ly) ; print("Ny = "+Ny) ; print("dy = "+dy) dy = Ly/Ny dx = Lx/Nx dt = 0.375 t = 0 ; Just for initial time psi(t,:,:) = fft2db(-fft2df(Zeta(0,:,:)*K2)) ; ##### Types of time differencing ##### ; time_diff = 1 ;Forward differencing ;time_diff = 2 ; Leap Frog ;time_diff = 3 ; Runge Kutta ; ##### Declare new arrays ##### fwd = new((/Tmax, Nx, Ny/), "float") dpsi_dy = new((/Tmax, Nx, Ny/), "float") dpsi_dx = new((/Tmax, Nx, Ny/), "float") u = new((/Tmax, Nx, Ny/), "float") v = new((/Tmax, Nx, Ny/), "float") dzeta_dy = new((/Tmax, Nx, Ny/), "float") dzeta_dx = new((/Tmax, Nx, Ny/), "float") ; if time_diff.eq.1 ; do t = 1, Tmax/dt print("t: "+t) ; fwd(t,:,:) = fft2db(-fft2df(Zeta(t,:,:))*K2) ; printVarSummary(fwd) ;psi(t,:,:) = fft2db(fwd(t,:,:)); ; end do ; end if ;Need to use Centered Finite Spatial Differencing (CD) to ;get u and v ; ##### For u=-d(psi)/dy using CD ##### do i = 0, Nsteps-1 print("i: "+i) if i.eq.0 dpsi_dy(t, i, :) = (psi(t,1,:) - psi(t, Ny-1,:))/(2*dy) else if i.gt.0.and.i.lt.Nsteps-1 dpsi_dy(t,i,:) = ((psi(t,i+1,:) - psi(t,i-1,:))/(2*dy)) else if i.eq.Nsteps-1 dpsi_dy(t,i,:) = (psi(t,0,:) - psi(t,Ny-2,:))/(2*dy) end if end if end if end do u(t,:,:) = -dpsi_dy(t,:,:) ; ##### For v = d(psi)/dx using CD ##### do j = 0, Nsteps-1 print("j: "+j) if j.eq.0 dpsi_dx(t,:,j) = (psi(t,:,1) - psi(t,:,Nx-1))/(2*dx) else if j.gt.0.and.j.lt.Nsteps-1 dpsi_dx(t,:,j) = (psi(t,:,j+1) - psi(t,:,j-1))/(2*dx) else if j.eq.Nsteps-1 dpsi_dx(t,:,j) = (psi(t,:,0) - psi(t,:,Nx-2))/(2*dx) end if end if end if end do v(t,:,:) = dpsi_dx(t,:,:) ; ##### For d(Zeta)/dY using CD ##### do i = 0, Nsteps-1 if i.eq.0 dzeta_dy(t,i,:) = (Zeta(t,1,:) - Zeta(t,Ny-1,:))/(2*dy) else if i.lt.0.and.i.gt.Nsteps-1 dzeta_dy(t,i,:) = (Zeta(t,i+1,:) - Zeta(t,i-1,:))/(2*dy) else if i.eq.Nsteps-1 dzeta_dy(t,i,:) = (Zeta(t,0,:) - Zeta(t,Ny-2,:))/(2*dy) end if end if end if end do ; ##### For d(Zeta)/dx using CD ##### do j = 0, Nsteps-1 if j.eq.0 dzeta_dx(t,:,j) = (Zeta(t,:,1) - Zeta(t,:,Nx-1))/(2*dx) else if j.lt.0.and.j.gt.Nsteps-1 dzeta_dx(t,:,j) = (Zeta(t,:,j+1) - Zeta(t,:, j-1))/(2*dx) else if j.eq.Nsteps-1 dzeta_dx(t,:,j) = (Zeta(t,:,0) - Zeta(t,:,Nx-2))/(2*dx) end if end if end if end do ;tst = fft2db(sf) ;printVarSummary(sf) ;wks = gsn_open_wks("x11","gyj") ;plot = gsn_csm_contour(wks,u,False) ;print(u) ;printVarSummary(v) ; ceof = fft2df(Zeta(0,:,:)) ; printVarSummary(ceof) ; xNew = fft2db(ceof) 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") ;gsn_define_colormap(wks,"WhViBlGrYeOrReWh") res = True res@cnFillOn = True res@cnMonoFillColor = False res@cnLineLabelsOn = False res@cnInfoLabelOn = False res@cnLinesOn = False ; res@mpGridSpacingF = 10. res@mpFillOn = False res@cnFillDrawOrder = "PreDraw" res@mpDataBaseVersion = "MediumRes" ; res@mpNationalLineThicknessF = 1.5 ; res@mpGeophysicalLineThicknessF = 2 ;For ShadeData res@cnLevelSelectionMode = "AutomaticLevels" res@tmXBLabelFontHeightF = fontHeightF res@tmYLLabelFontHeightF = fontHeightF res@gsnLeftStringFontHeightF = fontHeightF res@gsnRightStringFontHeightF = fontHeightF res@lbTitleFontHeightF = fontHeightF res@lbLabelFontHeightF = fontHeightF res@cnInfoLabelFontHeightF = fontHeightF res@lbBoxLinesOn = True res@lbLabelBarOn = True res@gsnLeftString = "Barotropic Model" res@gsnRightString = "Mike Ventrice" res@gsnSpreadColors = True res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = - 1 * 10^-3 res@cnMaxLevelValF = 1 * 10^-3 res@cnLevelSpacingF = 1 * 10^-4 ;plot = gsn_csm_contour( wks, Zeta(0,:,:), res ) plot = gsn_csm_contour( wks, psi(0,:,:), res ) ; Calculate Relative Vorticity ; vr = u ; dv = u ; av = u ; dav_dy = u ; uv2vrdvf(u,v,vr,dv) ; vr = vr*10^5 ; PI = acos(-1) ; f = 2 * 7.292e-5 * sin( lat * PI / 180 ) ; fConform = conform(vr, f, 1) ; av = vr + fConform ; y = lat * 111000 ; dav_dy = center_finite_diff_n(av, y, True, 0, 1) if( ( plotType.eq."png" ).or.( plotType.eq."gif" ) ) then system( "convert -trim -density 144 " \\ + plotName + ".eps " + plotName + "." + plotType ) system( "rm -f " + plotName + ".eps" ) end if print( "Plot is finished") print( "Thank you, come again." ) system( "date" ) end