; ; climo_1.ncl ; ; Concepts illustrated: ; - Calculating seasonal means ; - Creating a difference plot ; - Selecting specific decades of data ; - Calculating decadal means and standard deviations ; - Calculating probabilities ; - Calculating differences between decades ; - Selectively shading between contour levels (using an old method) ; - Overlaying a stipple pattern to show area of interest ; - Changing the density of contour shaded patterns ; - Drawing the zero contour line thicker ; - Copying coordinate arrays from one variable to another ; ;******************************************************** ; (1) Read NCEP/NCAR Reanalysis SLP ; (2) Use "runave" to compute seasonal [3-month] means ; (3) Use the "ind" function to select the 70s and 90s decades ; (4) Use clmMonLLT and stdMonLLT to compute the decadal ; means and standard deviations [contributed.ncl] ; (5) Use "ttest" to compute the probabilities ; (6) Use copy_VarCoords [contributed.ncl] to copy coordinate info ; (7) Calculate the differences between the decades ; (8) plot: ZeroNegDashLineContour and ShadeLtContour ; are in shea_util.ncl. ;*********************************************************** ; ; These files are loaded by default in NCL V6.2.0 and newer ; 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" ;A function for averaging june july august september undef ("season_from_month") function season_from_month (xMon:numeric, seastart:numeric, seaend:numeric, season:string) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; This function written by Gibies George is adapted from ; D Shea's month_to_seasonN() and month_to_season12() functions, ; with some slight difference in script, but enterly different usage ; which enables objective definision of season. ; ; Compute the seasonal average for user specified seasons. ; ; xMon(time) or xMon(time,lat,lon) or xMon(time,lev,lat,lon) ; ; The input "x" are assumed to contain monthly mean data ; The size of "time" MUST be divisible by 12. ; Also, it is assumed the "Jan" is the 1st month. ; ; first DJF season is a 2-month average (DJF=JF) ; ; USAGE : xSea = season_from_month (xMon, 6, 9, JJAS) ; gives "JJAS : Asian Summer Monsoon Season" ; : xSea = season_from_month (xMon, 10, 12, OND) ; gives "OND : Asian Post-monsoon Season" ; : xSea = season_from_month (xMon, 4, 5, pre-monsoon) ; gives "April-May : Asian Pre-monsoon Season" ; ; RESULT xSea(time/12,lat,lon) or xSea(time/12,lev,lat,lon) ; ; The above would return: ; xSea(time/12,lat,lon) or xSea(time/12,lev,lat,lon) ; ; NOTE: the "time" dimension may have to be altered to the user's desires. ; it may correspond to those associated with the 1st month. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; local dimx, rank, nmos, ntim, nyrs, xSea12, nlat, mlon \ , xSeaN, dName, cv begin isdef=isdefined((/"dimsize","dims","rank"/)) if(isdef(0)) delete(dimsize) end if if(isdef(1)) delete(dims) end if if(isdef(2)) delete(rank) end if delete (isdef) dimsize=dimsizes(xMon) dims=getvardims(xMon) rank = dimsizes(dimsize) if (rank.eq.2 .or. rank.ge.5) then print ("contributed: season_from_month: rank="+rank) print ("----- rank currently not handled -----") end if nmos = 12 ntim = dimsize(0) ;print ("contributed: season_from_month: ntim="+ntim+" nmos="+nmos) modCheck ("season_from_month", ntim, nmos) nyrs = ntim/nmos if (seaend.lt.seastart) then seaend = 12 + seaend end if sealen = seaend - seastart ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; xSea12 = month_to_season12 (xMon) ; compute the 12 seasons in month_to_seasonN() is here replaced ; by the following script adapted from month_to_season12(), ; with an additional the variable "sealen", ; which stands for the length of the season ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if (rank.eq.1) then ; (time) xSea = xMon ; transfer meta and reorder xSea = runave (xSea ,sealen, 0 ) ; overwrite with seasonal means xSea(0) = (xMon(0) + xMon(1) )*0.5 xSea(ntim-1) = (xMon(ntim-2) + xMon(ntim-1) )*0.5 xSea@long_name = season+" means: "+getLongName(xMon) xSea@season = season xSea12 = xSea($dims(0)$|:) end if if (rank.eq.3) then ; (time,lat,lon) ;xSea = xMon(lat|:,lon|:,time|:) ; transfer meta and reorder xSea = xMon($dims(1)$|:,$dims(2)$|:,$dims(0)$|:) ; transfer meta and reorder xSea = runave (xSea ,sealen, 0 ) ; overwrite with seasonal means xSea(:,:,0) = (xMon(0,:,:) + xMon(1,:,:) )*0.5 xSea(:,:,ntim-1) = (xMon(ntim-2,:,:) + xMon(ntim-1,:,:) )*0.5 xSea@long_name = season+" means: "+getLongName(xMon) xSea@season = season xSea12 = xSea($dims(0)$|:,$dims(1)$|:,$dims(2)$|:) end if if (rank.eq.4) then ; (time,lev,lat,lon) ;xSea = xMon(lev|:,lat|:,lon|:,time|:) ; transfer meta and reorder xSea = xMon($dims(1)$|:,$dims(2)$|:,$dims(3)$|:,$dims(0)$|:) xSea = runave (xMon ,sealen, 0 ) ; overwrite with seasonal means xSea(:,:,:,0) = (xMon(0,:,:,:) + xMon(1,:,:,:) )*0.5 xSea(:,:,:,ntim-1) = (xMon(ntim-2,:,:,:) + xMon(ntim-1,:,:,:) )*0.5 xSea@long_name = season+" means: "+getLongName(xMon) xSea@season = season xSea12 = xSea($dims(0)$|:,$dims(1)$|:,$dims(2)$|:,$dims(3)$|:) end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if (rank.ge.3) then nlat = dimsize(rank-2) mlon = dimsize(rank-1) end if if (rank.eq.1) then ; (time,lat,lon) xSeaN = new ( nyrs, typeof(xSea12), getFillValue(xMon)) ;printVarSummary(xSea12) xSeaN(:) = (/ xSea12(seastart:ntim-1:nmos) /) end if if (rank.eq.3) then ; (time,lat,lon) xSeaN = new ( (/nyrs,dimsize(1),dimsize(2)/), typeof(xSea), \ getFillValue(xMon)) ;printVarSummary(xSea12) xSeaN(:,:,:) = (/ xSea12(seastart:ntim-1:nmos,:,:) /) end if if (rank.eq.4) then ; (time,lat,lon) xSeaN = new ( (/nyrs,dimsize(1),dimsize(2),dimsize(3)/), typeof(xSea), \ getFillValue(xMon)) ;printVarSummary(xSea12) xSeaN(:,:,:,:) = (/ xSea12(seastart:ntim-1:nmos,:,:,:) /) end if if (rank.eq.5) then ; (time,lat,lon) xSeaN = new ( (/nyrs,dimsize(1),dimsize(2),dimsize(3),dimsize(4)/), typeof(xSea), \ getFillValue(xMon)) ;printVarSummary(xSea12) xSeaN(:,:,:,:,:) = (/ xSea12(seastart:ntim-1:nmos,:,:,:,:) /) end if copy_VarAtts (xMon, xSeaN) xSeaN@season = season if (isatt(xMon,"long_name") .or. isatt(xMon,"description") .or. \ isatt(xMon,"standard_name") ) then xSeaN@long_name = season+" Means: "+getLongName(xMon) end if ; copy dimension stuff dimName = xSea12!0 xSeaN!0 = dimName if(iscoord(xSea,dimName)) then cv = xSea12&$dimName$(seastart:ntim-1:nmos) xSeaN&$dimName$ = cv ; possibly override if (isatt(cv,"units") .and. \ (cv@units.eq."YYYYMM" .or. cv@units.eq."YYMM")) then cv = cv/100 cv@units = "YYYY" xSeaN&$dimName$ = cv end if if (isatt(cv,"units") .and. cv@units.eq."YYYYMMDD") then cv = cv/10000 cv@units = "YYYY" xSeaN&$dimName$ = cv end if end if if (rank.gt.1) then do i=1,rank-1 ; copy spatial coords dimName = xSea12!i xSeaN!(i) = dimName if(iscoord(xSea,dimName)) then xSeaN&$dimName$ = xSea12&$dimName$ end if end do end if ;printVarSummary(xSeaN) return (xSeaN) end ;============================================================================ ; This file still has to be loaded manually load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" ; ------------------------------------------------------- ;================================================================= begin ; ============================================================== ; User defined parameters that specify region of globe and ; ============================================================== ; ============================================================== ; Open the file: Read only the user specified period ; ============================================================== f = addfile ("hadsst_1950-2019.nc", "r") TIME = f->time date = cd_calendar(TIME,1) ; entire file sst = f->sst printVarSummary(sst) ; variable overview lat = f->latitude lon = f->longitude ; ============================================================== ; compute desired global seasonal mean: month_to_season (contributed.ncl) ; ============================================================== SST = season_from_month (sst, 10, 12, "DJF") nyrs = dimsizes(SST&time) printVarSummary(SST) ;******************************** ; indices for special dates ;********************************* ind5001 = ind(date.eq.195001) ; Jan 1950 ind8012 = ind(date.eq.198012) ; Dec 1980 ind8101 = ind(date.eq.198101) ; Jan 1981 ind1912 = ind(date.eq.201912) ; Dec 2019 ;******************************** ; compute [2-mon] climatologies ; over specified periods ;******************************** slpAve7079 = clmMonLLT ( SST(:,:,{ind5001:ind8012}) ) slpStd7079 = stdMonLLT ( SST(:,:,{ind5001:ind8012}) ) slpAve9099 = clmMonLLT ( SST(:,:,{ind8101:ind1912}) ) slpStd9099 = stdMonLLT ( SST(:,:,{ind8101:ind1912}) ) ;******************************** ; compute probabilities for means difference ;******************************** prob = ttest(slpAve9099, slpStd9099^2, 10 \ ,slpAve7079, slpStd7079^2, 10 ,True, False ) copy_VarCoords (slpAve9099, prob) prob@long_name = "Probability: difference between means" ; compute differences difAve = slpAve9099-slpAve7079 copy_VarCoords (slpAve9099, difAve) difAve@long_name = "SST-Difference" difAve@units = "mb" ;******************************** ; create plot ;******************************** nmo = 0 ; for demo, only plot Dec-Feb wks = gsn_open_wks ("png", "climo" ) ; send graphics to PNG file res = True ; plot mods desired res@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels res@cnMinLevelValF = -2. ; set min contour level res@cnMaxLevelValF = 2. ; set max contour level res@cnLevelSpacingF = 1. ; set contour spacing res@gsnDraw = False ; Do not draw plot res@gsnFrame = False ; Do not advance frome res@tiMainString = "SST Difference: " res@gsnCenterString = "5% stippled" res@gsnLeftString = "DJF" plot = gsn_csm_contour_map(wks,difAve(:,:,nmo), res) plot = ZeroNegDashLineContour (plot) res2 = True ; res2 probability plots res2@gsnDraw = False ; Do not draw plot res2@gsnFrame = False ; Do not advance frome res2@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels res2@cnMinLevelValF = 0.00 ; set min contour level res2@cnMaxLevelValF = 1.05 ; set max contour level res2@cnLevelSpacingF = 0.05 ; set contour spacing res2@cnInfoLabelOn = False res2@cnLinesOn = False ; do not draw contour lines res2@cnLineLabelsOn = False ; do not draw contour labels res2@cnFillScaleF = 0.6 ; add extra density delete(prob@long_name) ; add cyclic point plot2 = gsn_csm_contour(wks,gsn_add_cyclic_point(prob(:,:,nmo)), res2) opt = True ; set up parameters for pattern fill opt@gsnShadeFillType = "pattern" ; specify pattern fill opt@gsnShadeLow = 17 ; stipple pattern plot2 = gsn_contour_shade(plot2, 0.071, 300, opt) ; stipple all areas < 0.07 contour overlay (plot, plot2) draw (plot) frame(wks) end