load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" begin ;------------------------------------------------------------- ; pdo.create.ens.ncl -- script used to generate an ensenble PDO and related time series from different data sets ; C. Smith, NOAA/ESRL/PSD, 12/2018 ; PDO is defined as the leading EOF of monthly SSTA in the North Pacific, between 20N-70N ; SSTA = SST anomaly, where first the SST monthly climatology is removed and then the global mean is removed ; Output available at https://psl.noaa.gov/pdo/, where all PDOs are defined from 1920-2014 data ; NCL is available for download from https://www.ncl.ucar.edu/Download/ ;------------------------------------------------------------- ; MAIN CODE ;------------------------------------------------------------- ;------------------------------------------------------------- ; location of SST dataset(s) and land mask ; set to 0 if creating eofs ;------------------------------------------------------------- iupdate = 0 paths_ts = (/"/Datasets/COBE2/sst.mon.mean.nc","/Datasets.private/hadley/sst.mnmean.nc","/Datasets/noaa.ersst.v5/sst.mnmean.nc"/) path_landmask = (/"$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc"/) ; dataset names and associated plot labels dset_name=(/"COBE2-SST","HADISST1.1","NOAA ERSSTV5"/) top_lab = "Ensemble SST" ps_nameens = "sstens" ;------------------------------------------------------------- ; Time period for PDO EOF calculation (for PSD web page, use 192001 to 201412) ;------------------------------------------------------------- had = addfile(paths_ts(1),"r") cobe = addfile(paths_ts(0),"r") timecobe = cobe->time ctime = dimsizes(timecobe) timehad = had->time ntim = dimsizes(timehad) datehad = cd_calendar(timehad, 0) datecobe = cd_calendar(timecobe, 0) yearhad = tointeger(datehad(ntim-1,0)) monthhad = tointeger(datehad(ntim-1,1)) cmhad=sprinti("%0.2i", monthhad) yearcobe = tointeger(datecobe(ctime-1,0)) monthcobe = tointeger(datecobe(ctime-1,1)) cmcobe=sprinti("%0.2i", monthcobe) ymStrt = 192001 ; year-month start ymLast = 201412 ; year-month last ymStrto = 187001 ; year-month start earliest dataset cdate = yearcobe + cmcobe cdate = yearhad + cmhad ymLasto = cdate ; year-month last print(ymLasto +" ymLasto") year1 = 1870 year2 = yearhad ;------------------------------------------------------------- ; initalize variables ;------------------------------------------------------------- TimeDate = systemfunc("date") DATE = str_split(TimeDate, " ") UTC = DATE(1)+" "+ DATE(2)+" "+ DATE(5) nsim = dimsizes(paths_ts) print(nsim+" nsim") rad = 4*atan(1.0)/180 ;------------------------------------------------------------- ; common grid to be used for EOF calculation for all datasets ;------------------------------------------------------------- zlat = fspan(90.,-90.,91) zlat!0 = "lat" nlat = 91 zlon = fspan(0.,358.,180) zlon!0 = "lon" nlon = 180 nlatEOF = 26 latmin = 20. latmax = 70. nlonEOF = 76 lonmin = 110. lonmax = 260. description = new( 10, "string") file1_pdo_mon = new((/1,nlat,nlon/),"float") ;------------------------------------------------------------- ; loop through sst datasets ;------------------------------------------------------------- psd_txt_out = "updates/pdo.timeseries."+ps_nameens+".data" psd_csv_out = "updates/pdo.timeseries."+ps_nameens+".csv" eof_out = "pdo.eof."+ps_nameens+".nc" ts_out = "updates/pdo.timeseries."+ps_nameens+".nc" do nameind = 0,nsim-1 ; read this sst dataset sstdatafile = addfile(paths_ts(nameind),"r") sstall_in = sstdatafile->sst ; all of data TIME = sstdatafile->time sstall_in@missing_value=sstall_in@_FillValue ; ; regrid sst onto zlon/zlat (2x2) grid ; sstall = area_hi2lores_Wrap (sstall_in&lon,sstall_in&lat,sstall_in, True, 1, zlon,zlat, False) sstall!1 = "lat" sstall!2 = "lon" ; ;------------------------------------------------------------- ; now apply land mask to sst ;------------------------------------------------------------- d = addfile(path_landmask,"r") basemap = d->LSMASK(:,:) sstall = where(sstall.le.-1.8,-1.8,sstall) lsm = landsea_mask(basemap,sstall&lat,sstall&lon) sstall = mask(sstall,conform(sstall,lsm,(/1,2/)).ge.1,False) delete(lsm) delete([/basemap/]) ;------------------------------------------------------------- ; ; get time variables, index for start and end for EOF calculation ; ;------------------------------------------------------------- YYYYMM = cd_calendar(TIME, -1) ; convert YYYYMM1 = cd_calendar(TIME, 0) ; convert numtimes =dimsizes(TIME) iStrt = ind(ymStrt.eq.YYYYMM) ; index of start time iLast = ind(ymLast.eq.YYYYMM) ; last time iStrto = ind(ymStrto.eq.YYYYMM) ; index of start time iLasto = ind(ymLasto.eq.YYYYMM) ; last time year_all = tointeger(YYYYMM1(:,0)) date_start_yr = year_all(0) date_end_yr = year_all(numtimes-1) numtimeso = (iLasto -iStrto)+1 ;------------------------------------------------------------- ; subset of sst, over the EOF time period only, into sst_baseperiod ;------------------------------------------------------------- print(iStrt+" "+iLast+" "+ymStrt+" "+ymLast) sst_baseperiod = sstall(iStrt:iLast,:,:) print(iStrto+" "+iLasto+" "+ymStrto+" "+ymLasto) sst_overlap = sstall(iStrto:iLasto,:,:) sst_baseperiod!1 = "lat" sst_baseperiod!2 = "lon" ;------------------------------------------------------------- ; Now get SSTA (SST with annual cycle removed and then global mean SST each month removed) ; ******** ; First, remove annual cycle, defined only over base period, from SSTs ;------------------------------------------------------------- sstClm = clmMonTLL( sst_baseperiod ) ; get climo sst_baseperiod = calcMonAnomTLL (sst_baseperiod,sstClm) ; anomalies ;------------------------------------------------------------- ; remove same annual cycle from all ssts ;------------------------------------------------------------- kmm = 0 lyy = 1870 lyy = 1870 do it = 0,numtimeso-1 sst_overlap(it,:,:) = sst_overlap(it,:,:) - sstClm(kmm,:,:) kmm = kmm + 1 if(kmm.eq.12)then kmm = 0 lyy = lyy+1 end if end do ;------------------------------------------------------------- ; Then, remove area-averaged global mean SST each month ;------------------------------------------------------------- SSTA = sst_overlap ; initialize SSTA array SSTA_baseperiod = sst_baseperiod ; initialize SSTA_baseperiod array coswgt = cos(rad*sst_baseperiod&lat) ; cos weights coswgt!0 = "lat" coswgt&lat = sst_baseperiod&lat do ff = 0,dimsizes(sst_baseperiod&time)-1 SSTA_baseperiod(ff,:,:) = (/ sst_baseperiod(ff,:,:) - wgt_areaave(sst_baseperiod(ff,{-60:70},:),coswgt({-60.:70.}),1.0,0) /) SSTA(ff,:,:) = (/ sst_overlap(ff,:,:) - wgt_areaave(sst_overlap(ff,{-60:70},:),coswgt({-60.:70.}),1.0,0) /) end do ;------------------------------------------------------------- ; NOW AVERAGE; SET INITIAL ARRAYS FOR WHEN ON SST DATASET ;------------------------------------------------------------- if(nameind.eq.0)then indval = dimsizes(SSTA) ii = indval(0) jj = indval(1) kk = indval(2) indvalb = dimsizes(SSTA_baseperiod) iib = indvalb(0) jjb = indvalb(1) kkb = indvalb(2) SSTA_avg = new((/nsim,ii,jj,kk/),typeof(SSTA)) SSTA_base_avg = new((/nsim,iib,jjb,kkb/),typeof(SSTA_baseperiod)) end if SSTA_avg(nameind,:,:,:)=SSTA SSTA_base_avg(nameind,:,:,:)=SSTA_baseperiod print(SSTA_avg(nameind,0,45,90)+"test") print(SSTA_avg(nameind,ii-1,45,90)+"test") delete(sstall) delete(sstall_in) delete(SSTA) delete(sst_baseperiod) delete(TIME) delete( YYYYMM) delete( YYYYMM1) delete(year_all) end do SSTA = dim_avg_n_Wrap(SSTA_avg,0) SSTA_baseperiod=dim_avg_n_Wrap(SSTA_base_avg,0) ;------------------------------------------------------------- ; Get EOFs/PCs from SSTA ;------------------------------------------------------------- neofs = 1 ; number of EOFs output (could do more if desired) SSTA_baseperiod_CW = SqrtCosWeight(SSTA_baseperiod) ; Sqrt Cos Weight by latitude the SST array if(iupdate.eq.0)then evecv = eofunc_Wrap(SSTA_baseperiod_CW({lat|latmin:latmax},{lon|lonmin:lonmax},time|:),neofs,75) evecv!0 = "index" ; since default variable name is evn ; ; OUTPUT EOF TO NETCDF FILE print("fluffy nerd "+eof_out) system("/bin/rm -f "+eof_out) ; rm any pre-existing file eof = addfile(eof_out,"c") eof->eof = evecv else print("fluffy "+eof_out) eof = addfile(eof_out,"r") evecv = eof->eof end if ;------------------------------------------------------------- ; get PCs (PDO time series over base period) ; NOW need to get time series. ; keep SSTA_baseperiod_CW what is this? area weighted baseperiod ;------------------------------------------------------------- pcts_baseperiod = eofunc_ts_Wrap(SSTA_baseperiod_CW({lat|latmin:latmax},{lon|lonmin:lonmax},time|:),evecv,False) ; compute monthly principal component timeseries pcvari = tofloat(sprintf("%4.0f", evecv@pcvar(0))) ; percent variance explained delete(SSTA_baseperiod_CW) ; ; GET REGRESSION OF BASE PERIOD SSTA ON PDO TIME SERIES, put into pdo variable for later plotting pc1_std = dim_standardize(pcts_baseperiod(0,:),0) ; standardize the first PC (PDO PC) ; print(pc1_std) pdo = SSTA_baseperiod(0,:,:) ; initialize pdo array pdo&lat@units = "degrees_north" pdo&lon@units = "degrees_east" pdo = pdo@_FillValue pdo = (/ regCoef_n(pc1_std,SSTA_baseperiod,0,0) /) ; regress monthly SSTs on standardized PC1 pdo@pcvar = pcvari ;------------------------------------------------------------- ; Now get EOF time series (PC) from entire data record ;------------------------------------------------------------- SSTA_PDOregion = SSTA({lat|latmin:latmax},{lon|lonmin:lonmax},time|:) ; subset data into just PDO region ; what is SSTA. SSTA_PDOregion@missing_value = SSTA@_FillValue SSTA_PDOregion@_FillValue = SSTA@_FillValue nmon = dimsizes(SSTA_PDOregion&time) ;------------------------------------------------------------- ; reformat SSTA and PDO eof to get projection for PC time series ;------------------------------------------------------------- eofvec=reshape(evecv(lat|:,lon|:,index|:),(/nlatEOF*nlonEOF,neofs/)) xvec=reshape(SSTA_PDOregion(:,:,:),(/nlatEOF*nlonEOF,nmon/)) pcts = new((/neofs,nmon/),"float") do itime=0,nmon-1 pcts(:,itime) = dim_sum_n( eofvec*conform(eofvec,xvec(:,itime),0), 0 ) ; (neof) gets PCs at this time end do pcts!1 = "time" pcts&time = SSTA_PDOregion&time pc1 = pcts(0,:) delete([/pcts,evecv/]) ;------------------------------------------------------------- ; PDO IS DEFINED AS POSITIVE WHEN ANOMALY IN EASTERN PACIFIC IS POSITIVE, BUT EOFS HAVE ARBITRARY SIGN ;------------------------------------------------------------- if (pdo({37},{200}).ge.0) then pdo = pdo*-1. pc1 = pc1*-1. end if ; 1850 1920 2014 std_pc=stddev(pc1(348:1476)) pc1=pc1/std_pc file1_pdo_mon(0,:,:) = pdo file1_pdo_mon_timeseries = pc1 command="/bin/rm -f "+ts_out system(command) ; remove any pre-existing file timeseries = addfile(ts_out,"c") ;------------------------------------------------------------- ; create global file attributes ;------------------------------------------------------------- fAtt = True ; assign file attributes fAtt@title = "PDO timeseries from "+top_lab fAtt@Conventions = "CF-1.2" fAtt@creation_date = systemfunc ("date") fAtt@References = "https://psl.noaa.gov/pdo/" fAtt@references = "https://psl.noaa.gov/pdo/" fAtt@comment = "calcuated using NCL code from NOAA PSL" fAtt@instituion = "NOAA PSL" fAtt@source = "calculated using NCL code from NOAA PSL" fAtt@history = "generated from NCL code from NOAA PSL" fAtt@datasets = "NOAA ERSST V5/JMA COBE SST/UKMet HadISST1.1" filedimdef(timeseries,"time",-1,True) fileattdef( timeseries, fAtt ) ; copy file attributes timeseries->pdo = pc1 ; pdo!0 = "time" ; assign named dimensions timeseries->pdo@units = "C" timeseries->pdo@long_name="Pacific Decadal Oscillation (PDO)" timeseries->pdo@dataset=top_lab timeseries->pdo@statistic="PC" timeseries->pdo@comment="SST interpolated to a 2x2 grid" timeseries->pdo@missing_value=pc1@_FillValue range=ymStrt +" to "+ymLast timeseries->pdo@eof_range=range timeseries->pdo@eof_method="1st eof weighted SST from 20-90N" ; now write time series out in various formats. mon1=1 mon2=12 dsets="" do iname=0,nsim-1 dsets=dsets+dset_name(iname)+" " end do filename= psd_txt_out filenamecsv=psd_csv_out description(0)="PDO from "+top_lab description(1)="https://psl.noaa.gov//pdo/" description(2)="Using EOF from 1920 to 2014 for N Pacific (see webpage)" description(3)="Date produced "+UTC description(4)="SST interpolated to a 2x2 grid " description(5)="Datasets used "+dsets description(6)="Datasets used: " system("/bin/rm -f "+filename) ; rm any pre-existing file system("/bin/rm -f "+filenamecsv) ; rm any pre-existing file print( mon1+" "+mon2+" "+year1+" "+year2) ierr=writets ( pc1, mon1,mon2,year1,year2,filename,filenamecsv,description) delete(SSTA_PDOregion) delete( eofvec) delete( xvec) delete(file1_pdo_mon_timeseries) delete([/pdo,pc1/]) ;================================================================================================ ; Graphics section, spatial plot only ;------------------------------------ wks_type = "ps" if (wks_type.eq."png") then wks_type@wkWidth = 1500 wks_type@wkHeight = 1500 end if wks = gsn_open_wks(wks_type,"updates/pdo.eof."+ps_nameens) gsn_define_colormap(wks,"ncl_default") res = True res@mpProjection = "WinkelTripel" res@mpGeophysicalLineColor = "gray42" res@mpPerimOn = False res@mpGridLatSpacingF = 90 ; change latitude line spacing res@mpGridLonSpacingF = 180. ; change longitude line spacing res@mpGridLineColor = "transparent" ; trick ncl into drawing perimeter res@mpGridAndLimbOn = True ; turn on lat/lon lines res@mpFillOn = False res@mpCenterLonF = 180. res@mpOutlineOn = True res@mpDataSetName = "Earth..4" res@gsnDraw = False res@gsnFrame = False res@vpYF = 0.95 res@vpHeightF = 0.3 res@vpXF = 0.2 res@vpWidthF = 0.6 ; res@cnFillMode = "RasterFill" res@cnLevelSelectionMode = "ExplicitLevels" res@cnLevels = fspan(-.65,.65,27) res@cnLineLabelsOn = False res@cnFillOn = True res@cnLinesOn = False res@lbLabelBarOn = False res@gsnLeftStringOrthogonalPosF = -0.05 res@gsnLeftStringParallelPosF = .005 res@gsnCenterStringOrthogonalPosF = 0.05 res@gsnRightStringOrthogonalPosF = -0.05 res@gsnRightStringParallelPosF = 0.96 res@gsnRightString = "" res@gsnLeftString = "" res@gsnLeftStringFontHeightF = 0.0175 res@gsnCenterStringFontHeightF = 0.0185 res@gsnRightStringFontHeightF = 0.0175 res@gsnLeftString = "" res@gsnRightString = "" res@gsnCenterString = "" ; res@gsnLeftString = "a)" res@gsnCenterString = "PDO "+top_lab plot = gsn_csm_contour_map(wks,file1_pdo_mon(0,:,:),res) fheight=.01 txres = True txres@txFontHeightF = fheight ; Set the font height label = "NOAA PSL" ; txFuncCode is ~ gsn_text_ndc(wks,label,0.90, .20 ,txres) panres = True panres@gsnMaximize = True panres@gsnPanelBottom = 0.05 panres@gsnPaperOrientation = "portrait" panres@gsnPanelLabelBar = True panres@gsnPanelYWhiteSpacePercent = 3.0 panres@pmLabelBarHeightF = 0.05 panres@pmLabelBarWidthF = 0.55 panres@lbTitleOn = False panres@lbBoxLineColor = "gray70" panres@lbLabelFontHeightF = 0.0125 panres@txFontHeightF = 0.016 panres@txString = "" gsn_panel(wks,plot,(/1,1/),panres) delete(wks) ; system("convert -density 144 -trim +repage -border 8 -bordercolor white -background white -flatten "+get_script_prefi; x_name()+"."+wks_type+" "+get_script_prefix_name()+".gif") end