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" load "taylor_diagram.ncl" load "taylor_metrics_table.ncl" load "./taylor_stats.ncl" load "./taylor_diagram_cam.ncl" begin ; ============================================================== ; global ; ============================================================== yrStrt = 1970 yrLast = 1999 startt = get_cpu_time() ;-- returns the CPU time used by NCL ; ============================================================== ; Open the file: Read only the user specified period first observations then model ; ============================================================== f = addfile("air.mon.mean.nc", "r") ; TIME = f->time YYYY = cd_calendar(TIME,-1)/100 ; entire file iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast) air = f->air(iYYYY,:,:) air = air(:,::-1,:) ; make latitudes monotonically increasing (NCL syntax) printVarSummary(air) ; (time, lat,lon) AIR = month_to_season (air, "JJA") aClm = dim_avg_n_Wrap( AIR, 0) printVarSummary(aClm) ; (lat,lon) ; ============================================================== ; compute desired global annual or monthly climatology (12,nlat,mlon) ;============================================================model1deg cam5 temp f3= addfile("T2M_cam5.nc", "r") ; TIME3 = f3->time YYYY3 = cd_calendar(TIME3,-1)/100 ; entire file iYYYY3 = ind(YYYY3.ge.yrStrt .and. YYYY3.le.yrLast) T4 = f3->TREFHT(iYYYY3,:,:) T4=T4-273.15 printVarSummary(T4) ; (time, lat,lon) ;T4 = lonFlip( T4 ) t4 = linint2_Wrap(T4&lon, T4&lat, T4, True, air&lon, air&lat, 0) printVarSummary(t4) T4w = month_to_season (t4, "JJA") t4Clm = dim_avg_n_Wrap( T4w, 0) printVarSummary(t4Clm) ; (lat,lon) ;==================================================================precip f4 = addfile ("precip.mon.total.v301.nc", "r") TIME4 = f4->time YYYY4 = cd_calendar(TIME4,-1)/100 ; entire file iYYYY4 = ind(YYYY4.ge.yrStrt .and. YYYY4.le.yrLast) pm = f4->precip(iYYYY4,:,:) pm = pm(:,::-1,:) ; make latitudes monotonically increasing (NCL syntax) printVarSummary(pm) ; (time, lat,lon) PM = month_to_season (pm, "JJA") pClm = dim_avg_n_Wrap( PM, 0) printVarSummary(pClm) ; (lat,lon) ;=================================================================== f7= addfile("TotalPPcam5.nc", "r") ; TIME3 = f7->time YYYY3 = cd_calendar(TIME3,-1)/100 ; entire file iYYYY3 = ind(YYYY3.ge.yrStrt .and. YYYY3.le.yrLast) P3 = f7->Total_PP(iYYYY3,:,:) printVarSummary(P3) ; (time, lat,lon) ;P3 = lonFlip( P3 ) printVarSummary(P3) ; (time, lat,lon) p3 = linint2_Wrap(P3&lon, P3&lat, P3, True, pm&lon, pm&lat, 0) printVarSummary(p3) p3w = month_to_season (p3, "JJA") p3Clm = dim_avg_n_Wrap( p3w, 0) printVarSummary(p3Clm) ; (lat,lon) ;===========================================================================================PSL f88= addfile("SLP_PHIS_Z3cam5.nc", "r") ; TIME461 = f88->time YYYY461 = cd_calendar(TIME461,-1)/100 ; entire file iYYYY461 = ind(YYYY461.ge.yrStrt .and. YYYY461.le.yrLast) slpp = f88->PSL(iYYYY461,:,:) printVarSummary(slpp) ; (time, lat,lon) slpp=slpp/100 printVarSummary(slpp) ; (time, lat,lon) ; pm1 = slpp(:,::-1,:) ; make latitudes monotonically increasing (NCL syntax) ;printVarSummary(pm1) ; (time, lat,lon) s3w = month_to_season (slpp, "JJA") s3Clm = dim_avg_n_Wrap( s3w, 0) printVarSummary(s3Clm) ; (lat,lon) printMinMax(s3Clm, 0) ;================================================================== f6 = addfile ("slp.mon.mean.nc", "r") TIME46 = f6->time YYYY46 = cd_calendar(TIME46,-1)/100 ; entire file iYYYY46 = ind(YYYY46.ge.yrStrt .and. YYYY46.le.yrLast) psl = f6->slp(iYYYY46,:,:); psl = psl(:,::-1,:) printVarSummary(psl) S3 = linint2_Wrap(psl&lon, psl&lat, psl, True, slpp&lon, slpp&lat, 0) printVarSummary(S3) PSL = month_to_season (S3, "JJA") pslClm = dim_avg_n_Wrap( PSL, 0) printVarSummary(pslClm) ; (lat,lon) printMinMax(pslClm, 0) ;============================================================================================QFLX f11 = addfile ("eva_cam5mm.nc", "r") TIME411 = f11->time YYYY411 = cd_calendar(TIME411,-1)/100 ; entire file iYYYY411 = ind(YYYY411.ge.yrStrt.and. YYYY411.le.yrLast) qfx = f11->evacam5(iYYYY411,:,:) printVarSummary(qfx) PM1 = month_to_season (qfx, "JJA") eClm = dim_avg_n_Wrap( PM1, 0) printVarSummary(eClm) ; (lat,lon) ;================================================================== f71= addfile("lhtfl.mon.mean.nc", "r") ; TIME41 = f71->time YYYY41 = cd_calendar(TIME41,-1)/100 ; entire file iYYYY41 = ind(YYYY41.ge.yrStrt .and. YYYY41.le.yrLast) lm = f71->lhtfl(iYYYY41,:,:) printVarSummary(lm) lm = lm(:,::-1,:) ; make latitudes monotonically increasing (NCL l31 = linint2_Wrap(lm&lon, lm&lat, lm, True, qfx&lon, qfx&lat, 0) printVarSummary(l31) ; (time, lat,lon) l3w = month_to_season (l31, "JJA") l3Clm = dim_avg_n_Wrap( l3w, 0) printVarSummary(l3Clm) ; (lat,lon) ;=================================================================== f75= addfile("SLP_PHIS_Z3cam5.nc", "r") ; TIME35 = f75->time YYYY35 = cd_calendar(TIME35,-1)/100 ; entire file iYYYY35 = ind(YYYY35.ge.yrStrt .and. YYYY35.le.yrLast) z3 = f75->Z3(iYYYY35,:,:,:) printVarSummary(z3) ; (time, lev,lat,lon) z31 = z3(:,19,:,:) ; make latitudes monotonically increasing (NCL printVarSummary(z31) ; (time, lat,lon) Z3w = month_to_season (z31, "DJF") Z3Clm = dim_avg_n_Wrap( Z3w, 0) printVarSummary(Z3Clm) ; (lat,lon) ;====================================================================================z3 f45 = addfile ("hgt.mon.mean.nc", "r") TIME45 = f45->time YYYY45 = cd_calendar(TIME45,-1)/100 ; entire file iYYYY45 = ind(YYYY45.ge.yrStrt .and. YYYY45.le.yrLast) HGT = f45->hgt(iYYYY45,12,:,:) printVarSummary(HGT) ; (time, lat,lon) HGT3 = HGT(:,::-1,:) ; make latitudes monotonically increasing (NCL HGT33 = linint2_Wrap(HGT3&lon, HGT3&lat, HGT3, True, z31&lon, z31&lat, 0) HG = month_to_season (HGT33, "JJA") HClm = dim_avg_n_Wrap( HG, 0) printVarSummary(HClm) ; (lat,lon) ;====================================================================================z3 ;Taylor diagram calculations Generally, the 'common grid' is the control grid. ;=================================================================================== case = (/ "CESM CAM5"/) nCase = dimsizes(case ) ; # of Cases [Cases] var = (/ "T2M(Wilmott)", "PRECIP(Wilmott)","SLP(NCEP)","QLFX(NCEP)", "Z3(NCEP)"/) nVar = dimsizes(var) ; # of Variables SEASONS = (/ "ANN" /) nSeason = dimsizes( SEASONS ) ;================================================================== ; create weights: sqrt(cos(lat)) [or sqrt(gw) ] CESM lat/lon grid ; ================================================================= rad = 4.*atan(1.)/180. lat1 = f->lat clat = sqrt( cos(rad*lat1) ) printVarSummary(clat) copy_VarCoords(lat1, clat) printVarSummary(clat) lat2 = f88->lat clat2 = sqrt( cos(rad*lat2) ) copy_VarCoords(lat2, clat2) printVarSummary(clat2) ;==============================taylorstats======================== stat_taylor = taylor_stats(t4Clm, aClm, clat, 0) stat_taylor2 = taylor_stats(p3Clm, pClm, clat, 0) stat_taylor3 = taylor_stats(s3Clm, pslClm, clat2, 0) stat_taylor4 = taylor_stats(eClm, l3Clm, clat2, 0) stat_taylor5 = taylor_stats(Z3Clm, HClm, clat2, 0) CA_cc = (/stat_taylor(0), stat_taylor2(0), stat_taylor3(0),stat_taylor3(0),stat_taylor5(0)/) ; pattern_cor CA_rat = (/stat_taylor(1), stat_taylor2(1), stat_taylor3(1), stat_taylor4(0),stat_taylor5(0)/) ; ratio of CA_bias = (/stat_taylor(2), stat_taylor2(2), stat_taylor3(2),stat_taylor4(0),stat_taylor5(0)/) ; bias of ;********************************** ; Put the ratios and pattern correlations into arrays for plotting ;********************************** ; arrays to be passed to taylor_diagram. It will calculate the x xnd y coordinates. ratio = new ((/nCase, nVar/), typeof(CA_cc) ) cc = new ((/nCase, nVar/), typeof(CA_cc) ) bias = new ((/nCase, nVar/), typeof(CA_cc) ) ratio(0,:) = CA_rat cc(0,:) = CA_cc bias(0,:) = CA_bias print(bias) print(cc) print(ratio) ;================================ ; PLOTS ;================================ res = True ; default taylor diagram res@Markers = (/16/) ; make all solid fill res@Colors = (/"red"/) res@varLabels = var res@caseLabels = case res@varLabelsYloc = 0.6 ; Move location of variable labels ;;[default 0.45] res@caseLabelsFontHeightF = 0.14 ; make slight larger [default=0.12 ] res@varLabelsFontHeightF = 0.011 ; make slight smaller [default=0.013] res@stnRad = (/0.5, 1.25 /) ; res@tiMainString = "Annual Normalized Taylor Diagram" wks = gsn_open_wks("png","TaylorCAM5jja") ; send graphics to PNG file plot = taylor_diagram_cam(wks,ratio,cc,bias,res) ;res@ccRays = (/ 0.4, 0.9 /) ; correllation rays ;res@ccRays_color = "LightGray" ; default is "black" ; res@centerDiffRMS = False ; RMS 'circles' ;res@centerDiffRMS_color = "LightGray" ; default is "black"/) end