; Program to calculate eddy kinetic energy ;*********************************************************** ; Load the libraries ;*********************************************************** 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/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" begin ;*********************************************** ;Load CAM3.1 output files and read in data ;*********************************************** tStrt = 122 tLast = 481 tens = (tLast - tStrt) + 1 nlat = 128 mlon = 256 nlev = 26 ; do y = 1,2 y = 1 if ( y .eq. 5) y = 6 end if if (y.eq. 42) then y = 43 end if if ( y .eq. 54) then y = 55 end if in = addfile("/archive/u1/uaf/basu/ENSO_exp/control/control_sims06-60/control/control"+y+".cam2.h0.1997-11-01-00000.nc","r") temp = in->T(tStrt:tLast,:,:,:) u = in->U(tStrt:tLast,:,:,:) v = in->V(tStrt:tLast,:,:,:) uu = u^2 printVarSummary(v) printMinMax(v,True) printVarSummary(u) printMinMax(u,True) ;****************************************************** ; Thickness of pressure columns ;****************************************************** hyam = in->hyam ; read from a file the mid-layer coef hybm = in->hybm ; read from a file hyai = in->hyai ; read from a file the mid-layer coef hybi = in->hybi ; read from a file ps = in->PS(tStrt:tLast,:,:) ; surface pressure [Pa] p0 = in->P0 dp = dpres_hybrid_ccm(ps,p0,hyai,hybi) dp@long_name = "pressure" dp@units = "Pa" printVarSummary(dp) printMinMax(dp, True) copy_VarCoords(v,dp) delete(ps) ;*************************************************** ; Converting the variables to pressure surfaces ;*************************************************** ps = in->PS(tStrt:tLast,:,:) ; surface pressure [Pa] p0mb = 0.01*in->P0 ; since ps is in Pa or [ f->P0] pnew = (/850,775,700/) ;desired output levels(hPa/mb) intyp = 1 kxtrp = True nlev = dimsizes(hyam) tbot = temp(:,nlev-1,:,:) phis = in->PHIS(tStrt:tLast,:,:) varflg = 0 u_p = vinth2p_ecmwf(u,hyam,hybm,pnew,ps,intyp,p0mb,1,kxtrp,varflg,tbot,phis) printVarSummary(u_p) printMinMax(u_p, True) v_p = vinth2p_ecmwf(v,hyam,hybm,pnew,ps,intyp,p0mb,1,kxtrp,varflg,tbot,phis) printVarSummary(v_p) printMinMax(v_p, True) uu_p = vinth2p_ecmwf(uu,hyam,hybm,pnew,ps,intyp,p0mb,1,kxtrp,varflg,tbot,phis) printVarSummary(u_p) printMinMax(u_p, True) ;******************************************************* ; Calculate the zonal mean of zonal and meridional wind ;******************************************************* v_z = dim_avg_n_Wrap(v_p,3) printVarSummary(v_z) printMinMax(v_z, True) u_z = dim_avg_n_Wrap(u_p,3) printVarSummary(u_z) printMinMax(u_z, True) uu_z = dim_avg_Wrap(uu_p) printVarSummary(uu_z) printMinMax(uu_z, True) ;****************************************************** ; Meridional and zonal wind deviation from zonal mean ;****************************************************** vh = new((/tens,dimsizes(pnew),nlat,mlon/),typeof(v)) do i = 0,255 vh(:,:,:,i) = v_p(:,:,:,i) - v_z(:,:,:) end do copy_VarCoords(v_p,vh) printVarSummary(vh) printMinMax(vh, True) uh = new((/tens,dimsizes(pnew),nlat,mlon/),typeof(v)) do i = 0,mlon-1 uh(:,:,:,i) = u_p(:,:,:,i) - u_z(:,:,:) end do copy_VarCoords(u_p,uh) printVarSummary(uh) printMinMax(uh, True) ;*********************************************** ; Calculating Eddy Kinetic Energy ;*********************************************** eke = new((/tens,dimsizes(pnew),nlat,mlon/),typeof(v)) g =9.81 ;m/s^2 eke = (uh)^2 + (vh)^2 copy_VarCoords(u_p,eke) printVarSummary(eke) printMinMax(eke,True) ubb = u (1,0,:,:) eke775 = eke(:,1,:,:) printVarSummary(eke775) printMinMax(eke775,True) ;******************************************************* ; Calculate band pass filter weights ;******************************************************* ;lanczos weights for 2-6 day bandpass using 6-hourly data ;number of timesteps in low pass cutoff -- in this case ten days low = 24 ;number of timesteps in high pass cutoff -- in this case 3 days high= 8. nwt = 359 ;A scalar indicating the total number of weights ;(must be an odd number; nwt >= 3). ;The more weights, the better the filter, but there is a ;greater loss of data. fca = 1./low ;A scalar indicating the cut-off frequency of the ideal ;high or low-pass filter: (0.0 < fca < 0.5) fcb = 1./high ;A scalar used only when a band-pass filter is desired. ;It is the second cut-off frequency (fca < fcb < 0.5). ihp = 2 ;A scalar indicating the low-pass filter: ;ihp = 0; high-pass ihp = 1; band-pass ihp = 2. nsigma = 1. ;A scalar indicating the power of the sigma factor (nsigma >= 0). ;Note: nsigma=1. is common. ; compute weights wts = filwgts_lancos (nwt, ihp, fca, fcb, nsigma) ; print(wts) ; print(wts@resp+" "+wts@freq) ;******************************************************** ; apply lanczos wts to daily data -> band pass filtering ;******************************************************** eke_f = wgt_runave_n_Wrap(eke775,wts,0,0) printVarSummary(eke_f) printMinMax(eke_f,True) print(eke775(120,100,240)) print(eke_f(120,100,240)) ;************************************************* ; Seasonal mean of EKE ;************************************************* eke_mon_avg = dim_avg_n_Wrap(eke_f,0) copy_VarCoords(ubb,eke_mon_avg) printVarSummary(eke_mon_avg) printMinMax(eke_mon_avg, True) ;******************************************************* ; Write monthly mean to ascii file ;******************************************************* ; fmtx = "256f10.4" ; opt_1 = True ; opt_1@fout = "eke_bp_775_djf_"+y ; write_matrix(eke_mon_avg,fmtx,opt_1) ; end do ;************************************************* ;Plot ;************************************************* res = True res@gsnMaximize = True res@gsnSpreadColors = True res@cnFillOn = True res@cnLinesOn = False res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = 0 res@cnMaxLevelValF = 300 res@cnLevelSpacingF = 2 res@lbLabelStride = 10 res@lbOrientation = "vertical" ;vertical label bar res@tmLabelAutoStride = True res@tmXTOn = False res@tmYROn = False res@mpGeophysicalLineThicknessF = 1.5 res@mpCenterLonF = 180 ;**************************************************** ; Plot ;**************************************************** wks = gsn_open_wks("x11","eke2") gsn_define_colormap(wks,"BlAqGrYeOrReVi200") plot = gsn_csm_contour_map_ce(wks,eke_mon_avg,res) end