; Zonal Wind ;*************************************************************** 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" ;*************************************************************** begin ;****************************************************************************** ; read in data ECHAM5 ;***************************************************************************** in = addfile ("...ECHAM5.nc","r") time_EMAC = in->time lat_EMAC = in->lat lon_EMAC = in->lon lev_EMAC = in->lev hyam = in->hyam normhyam = hyam/101325 hybm = in->hybm hyai = in->hyai hybi = in->hybi aps_tmp = in->aps_ave ; TIME, LAT, LON aps = aps_tmp(:,::-1,:) U_wind_tmp = in->um1_ave ; TIME, LEV, LAT, LON U_wind = U_wind_tmp(:,:,::-1,:) U_wind!0 = "time" U_wind!1 = "lev" U_wind!2 = "lat" U_wind!3 = "lon" U_wind&time = time_EMAC U_wind&lev = lev_EMAC U_wind&lat = lat_EMAC U_wind&lon = lon_EMAC U_wind&lat = U_wind&lat(::-1) aps!0 = "time" aps!1 = "lat" aps!2 = "lon" aps&time = time_EMAC aps&lat = lat_EMAC aps&lon = lon_EMAC aps&lat = aps&lat(::-1) ;******************************************************** ; interpolation to new grid ;******************************************************** newgrid = g2fsh_Wrap(U_wind,(/73,144/)) printVarSummary(newgrid) newgrid2 = g2fsh_Wrap(aps,(/73,144/)) printVarSummary(newgrid2) ;******************************************************** ; Average Time ;******************************************************** u_ave = dim_avg_n_Wrap(newgrid(:,:,::-1,:),0) printVarSummary(u_ave) aps_ave = dim_avg_n_Wrap(newgrid2(:,::-1,:),0) printVarSummary(aps_ave) ;******************************************************* ; variance ;******************************************************* Std = dim_variance_n_Wrap(newgrid(:,:,::-1,:),0) printVarSummary(Std) ;******************************************************* ; interpolation to pressure coordinates ;******************************************************* lev_p = (/1000.,925.,850.,700.,600.,500.,400.,300.,250.,200.,150.,100.,70./) ccrpress = vinth2p (u_ave, normhyam, hybm, lev_p, aps_ave, 1, 1013.25, 1, False) ccrpress!0 = "lev_p" ccrpress!1 = "lat" ccrpress!2 = "lon" ccrpress&lev_p@units="hPa" Std_mmd = vinth2p (Std, normhyam, hybm, lev_p, aps_ave, 1, 1013.25, 1, False) Std_mmd!0 = "lev_p" Std_mmd!1 = "lat" Std_mmd!2 = "lon" Std_mmd&lev_p@units="hPa" ;********************************************************* ; Average Longitude ;********************************************************* clim_zon_wind = dim_avg_Wrap(ccrpress(lev_p|:,lat|:,lon|:)) ; LEV, LAT clim_zon_wind!0="lev" clim_zon_wind!1="lat" Std_ave = dim_avg_Wrap(Std_mmd(lev_p|:,lat|:,lon|:)) ; LEV, LAT Std_ave!0="lev" Std_ave!1="lat" ;******************************************************************************************* ; read in data NCEP ;******************************************************************************************* setfileoption("nc","MissingToFillValue",False) in2 = addfile("...uwnd.mon.mean.nc","r") u2 = in2->uwnd(:,:,:,:) time_OBS = in2->time lat_OBS = in2->lat lon_OBS = in2->lon lev_OBS = in2->level u2!0 = "time" u2!1 = "lev" u2!2 = "lat" u2!3 = "lon" u2&time = time_OBS u2&lev = lev_OBS u2&lat = lat_OBS u2&lon = lon_OBS ;******************************************************************* ; Average Time ;******************************************************************* u2_ave = dim_avg_n_Wrap(u2(361:720,:,:,:),0) printVarSummary(u2_ave) ;******************************************************************* ; Variance ;******************************************************************* Std2 = dim_variance_n_Wrap(u2(361:720,:,:,:),0) printVarSummary(Std2) ;********************************************************************* ; new pressure levels ;********************************************************************* ; desired pressure levels: pnew = (/1000.,925.,850.,700.,600.,500.,400.,300.,250.,200.,150.,100.,70./) ; old pressure levels lev = (/1000.,925.,850.,700.,600.,500.,400.,300.,250.,200.,150.,100.,70.,50.,30.,20.,10./) lev@units = "mb" ; units attribute required by gsn_csm_pres_hgt ;*************************************************************************** ; interpolate von old pressure levels to new pressure levels ;*************************************************************************** clim_zon_wind2 = int2p_n_Wrap(lev,u2_ave(lat|:,lon|:,lev|:),pnew,2,2) clim_zon_wind2!0="lat" clim_zon_wind2!1="lon" clim_zon_wind2!2="pnew" clim_zon_wind2&lat = lat_OBS clim_zon_wind2&lon = lon_OBS clim_zon_wind2&pnew@units = "hPa" result = dim_avg_Wrap(clim_zon_wind2(pnew|:,lat|:,lon|:)) Std2_mmd = int2p_n_Wrap(lev,Std2(lat|:,lon|:,lev|:),pnew,2,2) Std2_mmd!0="lat" Std2_mmd!1="lon" Std2_mmd!2="pnew" Std2_mmd&lat = lat_OBS Std2_mmd&lon = lon_OBS Std2_mmd&pnew@units = "hPa" Std2_ave = dim_avg_Wrap(Std2_mmd(pnew|:,lat|:,lon|:)) ;************************************************************ ; Difference ;************************************************************ Delta = clim_zon_wind-result Delta!0 = "pnew" Delta!1 = "lat" Delta&pnew = pnew Delta&lat = lat_OBS Delta&pnew@units = "hPa" if(any(ismissing(Delta))) then print("Your data contains some missing values. Beware.") end if if (any(isnan_ieee(Delta))) then value = -1.e34 replace_ieeenan (Delta, value, 0) Delta@_FillValue = value end if printVarSummary(Delta) ;******************************************************************* ; T TEST ;******************************************************************* prob = 100.*(1. - ttest(result, Std2_ave, 30, clim_zon_wind, Std_ave, 30, True, False)) prob!0 = "lev_p" prob!1 = "lat" prob&lev_p = lev_p prob&lat = lat_OBS printVarSummary(prob) copy_VarCoords (Std2_ave, prob) ;********************************************************************************************* ; Plot ;********************************************************************************************* wks = gsn_open_wks("pdf","zonal_wind_0-OBS") gsn_define_colormap(wks,"BlWhRe") res = True res@cnLevelSelectionMode = "AutomaticLevels" res@gsnSpreadColors = True res@cnFillOn = True res@cnLineLabelsOn = True res@cnLineLabelBackgroundColor = -1 res@tiXAxisString = "Geographische Breite" res@tiYAxisString = "Druck in hPa" res@tiXAxisFontHeightF = 0.018 res@tiYAxisFontHeightF = 0.018 res@tmYRMode = "Automatic" res@lbOrientation = "vertical" res@lbLabelAutoStride = True ccrpresplot = gsn_csm_pres_hgt(wks,Delta(pnew|:,{lat|-90:90}),res) end