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"

 begin
; ==============================================================
; User defined parameters that specify region of globe and
; ==============================================================
latS   =  29.
latN   =  60.
lonL   = -15.
lonR   =  66.

yrStrt = 1900
yrLast = 1999

; ==============================================================
; Open the file: Read only the user specified period first observations then model
; ==============================================================
 season = "DJF"    ; choose Dec-Jan-Feb seasonal mean
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,:,:)
printVarSummary(air)                                ; (time, lat,lon)
lat = f->lat
lon= f->lon
;clat=cos(0.01745329*lat)           ; cos(lat)weight
; ==============================================================
; compute desired global seasonal mean: month_to_season (contributed.ncl)
; ==============================================================
air    = lonFlip( air )
AIR    = month_to_season (air, season)
printVarSummary(AIR)
nyrs   = dimsizes(AIR&time)
; =================================================================
; create weights:  sqrt(cos(lat))   [or sqrt(gw) ]
; =================================================================
  rad    = 4.*atan(1.)/180.
 lat = f->lat
   if (typeof(lat).eq."double") then
       clat = sqrt( cos(rad*tofloat(lat)) )
   else
       clat = sqrt( cos(rad*lat) )
   end if
   copy_VarCoords(lat, clat) ; contributed
   printVarSummary(clat) 
  ;clat = tofloat(f->lat) ; convert to float     
  ;clat   = sqrt( cos(rad*clat) )  ; gw for gaussian grid but clat is a double array
; =================================================================
; weight all observations 
; =================================================================
  wAIR   = AIR  ; type float
  wAIR   = AIR*conform(AIR, clat, 1);creates an array that is the same size and shape as SLP (wSLP) but is populated with entries of clat which is type double. The product of a float*double is a double. 
  printVarSummary(wAIR)                                  ; copy meta data
; =================================================================
; Reorder (lat,lon,time) the *weighted* input data
; Access the area of interest via coordinate subscripting
; =================================================================
airw     = AIR({lat|latS:latN},{lon|lonL:lonR},time|:)
  ;x      = wAIR(time|:,{lat|latS:latN},{lon|lonL:lonR})
printVarSummary(airw)
;airw = dtrend(airw(lat|:,lon|:,time|:),False)
printVarSummary(airw)

;============================================================model

f1= addfile("T2M_cam5.nc", "r") ;
TIME1   = f1->time
YYYY1   = cd_calendar(TIME1,-1)/100                 ; entire file
iYYYY1  = ind(YYYY1.ge.yrStrt .and. YYYY1.le.yrLast)
t2m    = f1->TREFHT(iYYYY1,:,:)
t2m=t2m-273.15
printVarSummary(t2m)                                ; (time, lat,lon)
;t2m    = lonFlip( t2m )
t2mw    = month_to_season (t2m, season)
printVarSummary(t2mw)
nyrs   = dimsizes(t2mw&time)
lat1 = f1->lat
 lon1= f1->lon
; =================================================================
; Reorder (lat,lon,time) the *weighted* input data
; Access the area of interest via coordinate subscripting
; =================================================================
;t2mw     = T2M({lat|latS:latN},{lon|lonL:lonR},time|:)
  ;x      = wT2M(time|:,{lat|latS:latN},{lon|lonL:lonR})
;printVarSummary(t2mw)
;t2mw = dtrend(t2mw(lat|:,lon|:,time|:),False)
;printVarSummary(t2mw)

;===============================================================
;==================================================================
;Taylor diagram calculations
;================================ intrapolation onto common grid (of observational data)
;ncl pattern_cor between different size arrays
;************************************************
; interpolate to new grid
;***********************************************
  ;newlat = fspan(-60.,60,24)
  ;newlon = fspan(0.,355.,72)

 newt2mw = linint2_Wrap(lon1,lat1,t2m,True,lon,lat,0)

  newt2mw!0   ="lat1"
  newt2mw!1   = "lon1"
  newt2mw&lat = lat
  newt2mw&lon = lon
;======================================================================centered Pattern correlation (coslat weighting has been done previously above)
re=escorc(airw,newt2mw)
 ;cor1 = dim_avg_n_Wrap(pattern_cor( t2mw, airw, clat, 0), 0);rc = pattern_cor(x, y,gw, 0)      ; gaussian weighting, centered
 mmd= (/cor1/)
 printVarSummary(mmd)
;================================Standard deviation ;================================

;pre0_Std = dim_avg_n_Wrap( dim_stddev_n_Wrap( t2mw, (/1,2/)), 0)
 ;std1 = dim_rmsd_Wrap(airw,t2mw, 0);computes rootmean square difference
;std1 = dim_rmsd_n(t2mw, airw, 0);
    std1 = dim_rmsd( t2mw(lat|:,lon|:,time|:), airw(lat|:,lon|:,time|:) )    ; ==> rmsdTime(nlat,nlon)

    ;rmsdTime = dim_rmsd_n( x, y, 0 )                                       ; ==> no reordering needed
;================================Standard deviations normalisation;================================
 ;std1 = t2mw_Std/airw_Std
;================================
; Cases [Model]
 case = (/ "precipitation" /)
 nCase = dimsizes(case )
; variables compared
var = (/ "CMAP"/)
 nVar = dimsizes(var)

 p_rat = (/std1/)
p_cc = (/cor1/)

;;================================ Put the ratios and pattern correlations into arrays for plotting
 nDataSets = 1 ; number of datasets
npts = dimsizes(p_rat)
; arrays to be passed to taylor plot 
 ratio = new ((/nCase, nVar/), typeof(p_rat))
cc = new ((/nCase, nVar/), typeof(p_cc) )

ratio(0,:) = p_rat
cc(0,:) = p_cc
;================================
; PLOTS
;============================================================
 res = True ;diagram mods desired
res@tiMainString = "precipitation" ; title
 res@Colors = (/"blue"/) ; ; marker colors
res@Markers = (/14/) ;,9,2,3,8/) ; marker styles
res@markerTxYOffset = 0.03 ;offset btwn marker & label
res@gsMarkerSizeF = 0.018 ;

 res@txFontHeightF = 0.018 ;
 res@stnRad = (/ 0.75, 1.25 /) ; additional standard radii
 res@ccRays = (/ 0.6, 0.9 /) ; correllation rays
 ;res@ccRays_color = "LightGray" ; default is "black"
 res@varLabels = var
 res@caseLabels = case ; affiche la liste des variables chargées
 res@varLabelsYloc = 1.5 ; Move location of variablelabels [pour taylor]
;res@varLabelsYloc = 1.65 ; Move location of variable labels [pour taylor2]
;res@varLabelsYloc = 2.04 ; Move location of variablelabels [pour taylor1]
res@caseLabelsFontHeightF = 0.10 ; make slight larger[default=0.12 ]
 res@varLabelsFontHeightF = 0.01 ; make slight smaller[default=0.013]
 res@centerDiffRMS = True ; RMS 'circles'
 res@centerDiffRMS_color = "LightGray" ; default is "black"
 wks = gsn_open_wks("eps", "pre_tailor_annual")
plot = taylor_diagram(wks, ratio, cc, res)
end