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 ; ****************** setting ************ nsvd = 2 year1 = 1979 year2 = 2015 nyear = year2-year1+1 year = ispan(year1,year2,1) ymStrt = year1*100 + 1 ymLast = year2*100 + 12 season = (/"DJF","MAM","JJA","SON"/) ns = dimsizes(season) ; area set for the SST data latS = -90. latN = -45. lonW = 0 lonE = 360 detrend = "" ; *************************** setting over ******** ; **************************** read data ********* ; read the NCEP SLP fi = addfile("/Users/linx/data/NCEP/reanalysis.derived/surface/slp.mon.mean.nc","r") time = fi->time yymm = ut_calendar(time,-1) nt1 = ind(yymm.eq.ymStrt) nt2 = ind(yymm.eq.ymLast) data = fi->slp(nt1:nt2,{latS:latN},{lonW:lonE}) slp = data*data@scale_factor + data@add_offset copy_VarMeta(data,slp) slp@_FillValue = slp@missing_value delete(time) delete(yymm) ; read sst fi = addfile("/Users/linx/data/Hadley/HadISST.2.2.0.0_sea_ice_concentration.nc","r") time = fi->time yymm = ut_calendar(time,-1) nt1 = ind(yymm.eq.ymStrt) nt2 = ind(yymm.eq.ymLast) sst = fi->sic(nt1:nt2,{-45:-80},{lonW:lonE}) printVarSummary(sst) ; ************************** read over ********************* ; ********* prepare the array before Loop *************** nlat = dimsizes(slp&lat) nlon = dimsizes(slp&lon) mlat = dimsizes(sst&latitude) mlon = dimsizes(sst&longitude) svd_slp_season = new((/ns,nsvd,nlat,nlon/),"float") svd_sst_season = new((/ns,nsvd,mlat,mlon/),"float") ts_slp_season = new((/ns,nsvd,nyear/),"float") ts_sst_season = new((/ns,nsvd,nyear/),"float") ; ************************* Loop ************************ do i = 0,0 ;ns-1 print(i + " " + season(i) + " ") ; 1. calculate the seasonal average and detrend slp_season = month_to_season(slp, season(i)) sst_season = month_to_season(sst, season(i)) ;printVarSummary(slp_season) if(detrend.eq."dtrend")then slp_ano = dtrend_msg(year, slp_season(lat|:,lon|:,time|:), False, False) sst_ano = dtrend_msg(year, sst_season(latitude|:,longitude|:,time|:), False, False) print(" " + "dtrend data") else slp_ano = slp_season(lat|:,lon|:,time|:) sst_ano = sst_season(latitude|:,longitude|:,time|:) print(" " + "no dtrend") end if ; 2. reshape slp yslp = reshape(slp_ano,(/nlat*nlon,nyear/)) ; 2. reshape sst ls = dim_avg_Wrap(sst_season(latitude|:,longitude|:,time|:)) nsst = num(ls.ne.sst_season@_FillValue) ; 14868 grids with non-missing values for all months print("total numbers of ocean grids in the chosen area : " + nsst) ls1d = ndtooned(ls) lsind = ind(ls1d.ne.sst_season@_FillValue.and.ls1d.ne.0.and.ls1d.ne.1) nsst = num(ls1d.ne.sst_season@_FillValue.and.ls1d.ne.0.and.ls1d.ne.1) print("after removing constant 0/1 value," + nsst + " grids left" ) ysst = new((/nsst,nyear/),"float") do j = 0,nyear-1 sst_1d = ndtooned(sst_ano(:,:,j)) ysst(:,j) = sst_1d( lsind ) end do ysst!0 = "grid" ysst!1 = "time" ; 3. calculate the SVD print(" " + num(ismissing(yslp))) print(" " + num(ismissing(ysst))) print(" " + num(yslp.eq.0)) print(" " + num(ysst.eq.0)) print(" " + num(ysst.eq.1)) ncols = nlat*nlon ncolz = nsst homlft = new((/nsvd,ncols/),float) hetlft = new((/nsvd,ncols/),float) homrgt = new((/nsvd,ncolz/),float) hetrgt = new((/nsvd,ncolz/),float) svd = svdstd(yslp, ysst, nsvd, homlft, hetlft, homrgt, hetrgt ) ;print( " svd done" ) print("svdstd: percent variance= " + svd) slp_ts = eofcov_ts( yslp, hetlft ) sst_ts = eofcov_ts( ysst, hetrgt ) ;printVarSummary(hetlft) ;printVarSummary(hetrgt) ;printVarSummary(svd) ; 4. reshape the output into original array svd_slp = reshape(hetlft,(/nsvd,nlat,nlon/)) svd_slp!0 = "nsvd" svd_slp!1 = "lat" svd_slp!2 = "lon" svd_slp&lat = slp&lat svd_slp&lon = slp&lon svd_tmp = reshape(hetrgt,(/nsvd,nsst/)) svd_sst = new((/nsvd,mlat,mlon/),"float") svd_sst = sst@_FillValue svd_sst@_FillValue = sst@_FillValue do k = 0,nsvd-1 sv1d = ndtooned(svd_sst(k,:,:)) sv1d(lsind) = svd_tmp(k,:) svd_sst(k,:,:) = onedtond(sv1d, (/mlat,mlon/)) end do svd_sst!0 = "nsvd" svd_sst!1 = "lat" svd_sst!2 = "lon" svd_sst&lat = sst&latitude svd_sst&lon = sst&longitude ; 5. arrage into new array svd_slp_season(i,:,:,:) = svd_slp svd_sst_season(i,:,:,:) = svd_sst ts_slp_season(i,:,:) = slp_ts ts_sst_season(i,:,:) = sst_ts svd_slp_season!0 = "season" svd_slp_season!1 = "nsvd" svd_slp_season!2 = "lat" svd_slp_season!3 = "lon" svd_slp_season&lat = slp&lat svd_slp_season&lon = slp&lon svd_sst_season!0 = "season" svd_sst_season!1 = "nsvd" svd_sst_season!2 = "lat" svd_sst_season!3 = "lon" svd_sst_season&lat = sst&latitude svd_sst_season&lon = sst&longitude ts_slp_season!0 = "season" ts_slp_season!1 = "nsvd" ts_slp_season!2 = "time" ts_slp_season&time = ispan(year1,year2,1) copy_VarCoords(ts_slp_season, ts_sst_season) delete(ls1d) delete(lsind) delete(ysst) delete(sst_1d) delete(sst_ano) delete(homrgt) delete(hetrgt) delete(svd_tmp) end do ; output netCDF if(detrend.eq."")then output_filename = "SVD_sic.vs.slp_seasonal_" + year1 + "-" + year2 else output_filename = "SVD_sic.vs.slp_seasonal.dtrend_" + year1 + "-" + year2 end if system("rm -f " + output_filename + ".nc") ; remove any pre-existing file ncdf = addfile(output_filename+".nc","c") ; open output netCDF file ;=================================================================== ; create global attributes of the file (optional) ;=================================================================== fAtt = True ; assign file attributes ; fAtt@title = "NSIDC Bootstrap sea ice concentration" ;fAtt@source_file = "original-file.nc" ;fAtt@Conventions = "None" fAtt@creation_date = systemfunc ("date") fileattdef( ncdf, fAtt ) ; copy file attributes ;=================================================================== ; make time an UNLIMITED dimension; recommended for most applications ;=================================================================== filedimdef(ncdf,"time",-1,True) ;=================================================================== ; output variables directly; NCL will call appropriate functions ; to write the meta data associated with each variable ;=================================================================== ncdf->svd_slp_season = svd_slp_season printVarSummary(svd_slp_season) printVarSummary(svd_sst_season) ncdf->svd_sst_season = svd_sst_season ncdf->ts_slp_season = ts_slp_season ncdf->ts_sst_season = ts_sst_season end