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" 
;****************************************************
 begin

; ==============================================================
; NAO region + Caspian region
; ==============================================================
  latS   =  30.
  latN   =  63. 
  lonL   = 30.
  lonR   =  63

  yrStrt = 1900
  yrLast = 1999

; ==============================================================
; 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)
lat = f->lat
lon= f->lon
  AIR=month_to_annual(air,1)                
      obst=AIR(:,0,0)  
         aClm = dim_avg_n_Wrap( AIR, 0)
   printVarSummary(aClm)                 ; (lat,lon)
   printMinMax(aClm, 0)               
; ==============================================================
; 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_annual(t4,1)                                           ;AIR = clmMonTLL (air)   
 modelt=T4w(:,0,0)
t4Clm = dim_avg_n_Wrap( T4w, 0)
   printVarSummary(t4Clm)                 ; (lat,lon)
   printMinMax(t4Clm, 0)  

  time = ispan(1900,1999,1)
  time@long_name = "Time"
printVarSummary(time)                                



  wks = gsn_open_wks("png","overlaytemp") ; send graphics to PNG file

; resources for "left" variable
   res = True
  colors = (/"red","blue"/)
  res@xyXStyle               = "Time"
  res@gsnMaximize      = True
  res@gsnPaperOrientation = "portrait"
  res@gsnDraw          = False
  res@gsnFrame         = False
  res@tiXAxisString          = "Year"
  res@xyLineThicknessF = 2.0
  res@xyLineColor      = colors(0)
  plot0 = gsn_csm_xy(wks,time,obst,res)
  res@xyLineColor     = colors(1)
  plot1 = gsn_csm_xy(wks,time,modelt,res)
  overlay(plot0,plot1)

lineres = True
  lineres@lgLineColors = (/"red","orange"/) ; line colors
  lineres@lgLineThicknesses = 3                        ; line thicknesses
  lineres@LineLengthPercent = 9.                         ; expressed as %, 0->100, length of line

  textres = True
  textres@lgLabels = (/"2mTemperature(Wilmott)","2mTemperature(CESM1degCAM5)"/)  ; legend labels (required)

  plot = simple_legend(wks,plot1,res,lineres,textres)

  draw(plot0)
  frame(wks)
end