;-----------------------------------------------------------
;    Dr. Md. Jalal Uddin
;    Postdoctoral Researcher, Jeju National University
;    PhD in Atmospheric Physics and Atmospheric Environment (NUIST)
;    M.Sc. in Applied Meteorology (NUIST)
;    B.Sc. in Disaster Management (PSTU)
;    Founder and Director of Research Society
;    Email:        dmjalal90@gmail.com, 20205103002@nuist.edu.cn    
;    Website:      https://researchsociety20.org/founder-and-director/ 
;    Researchgate: https://www.researchgate.net/profile/Md_Uddin125  
;-----------------------------------------------------------

; open data file

  f      = addfile ("sst.mon.anom.nc", "r")
 
  TIME   = f->time
  YYYY   = cd_calendar(TIME,-1)/100      ;cd_calendar (time, option) ;option = -1 or 1. The values returned will be in the format YYYYMM and will be type double for option = 1, and type integer for option = -1.        
  ;print(YYYY)
  
  
  ; define starting and ending year
  yrStrt = 2001
  yrLast = 2022
  
  iYYYY  = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)  ;Returns the indices where the input is True. A one-dimensional logical array.
  ;print(iYYYY)
  
  sst    = f->sst(iYYYY,:,:)
  printVarSummary(sst)  
  
  
  nyrs   = dimsizes(sst&time)
              
 ;----------------------------------------------
 ; EOF calculation

 ; create weights:  sqrt(cos(lat))   [or sqrt(gw) ]

  rad    = get_d2r("float")
  clat   = f->lat           
  clat   = sqrt(cos(rad*clat))        ; gw for gaussian grid

; weight all observations 

  wsst   = sst                                   ; copy meta data
  wsst   = sst*conform(sst, clat, 1)             ; conform (data,subset of sst,ndim)   https://www.ncl.ucar.edu/Document/Functions/Built-in/conform.shtml
  
  printVarSummary(wsst)
  
 
  ; Calculate EOFs (empirical orthogonal functions) of the SST over the Pacific Ocean.
  ; eofunc_n_Wrap:  https://www.ncl.ucar.edu/Document/Functions/Contributed/eofunc_n_Wrap.shtml
  ; eofunc_ts_n_Wrap: https://www.ncl.ucar.edu/Document/Functions/Contributed/eofunc_ts_n_Wrap.shtml


  neof         = 1        ; number of EOFs
  optEOF       = True       
  ;optEOF@jopt  = 0   ; The default is to use a covariance matrix 
  optEOF@jopt = 1   ; **only** if the correlation EOF is desired
  optETS      = False 

  
  eof    = eofunc_n_Wrap(wsst, neof, optEOF, 0)      ; eofunc_n_Wrap (data, neval,optEOF,dim )
  ;eof_ts = eofunc_ts_n_Wrap (wsst, eof, optETS, 0)   ; eofunc_ts_n_Wrap (data,evec,optETS,dim)

  printVarSummary(eof)                         ; examine EOF variables
  ;printVarSummary(eof_ts)


;------------------------------------------------
; creat plot
;------------------------------------------------

  pltType = "png"            
  pltDir  = "./"
  pltName = "eof_global_nino"
  pltPath = pltDir + pltName 
  
  pltType@wkWidth  = 2500
  pltType@wkHeight = 2500 

  wks  = gsn_open_wks(pltType,pltPath)

 
  res                      = True         
  res@gsnDraw              = False        ; don't draw yet
  res@gsnFrame             = False        ; don't advance frame yet
  res@gsnAddCyclic         = True        ; plotted dataa are not cyclic
  
  res@mpFillOn             = False        ; turn off map fill
  res@mpFillDrawOrder      = "PostDraw"          ; color of land
  res@mpLandFillColor      = "grey"              ; color of land
  res@mpGeophysicalLineThicknessF = 4
  res@mpGeophysicalLineColor      = "blue"
  res@mpNationalLineColor         = "blue"
  res@mpOutlineBoundarySets       = "AllBoundaries"
  res@pmTickMarkDisplayMode       = "Always"
   
     
  ; subset of region interest
  ;minLon = -270
  ;maxLon = -70
  ;minLat = -50
  ;maxLat = 50
  
  res@mpCenterLonF = 180.
  ;res@mpMinLatF    = minLat
  ;res@mpMaxLatF    = maxLat
  ;res@mpMinLonF    = minLon
  ;res@mpMaxLonF    = maxLon
  
  res@cnFillOn             = True         ; turn on color fill
  res@cnLinesOn            = False        ; True is default
  res@cnLineLabelsOn       = False            ; turn off line labels
  res@cnInfoLabelOn        = False            ; turn off info label
  
  res@cnFillMode           = "RasterFill"
  res@cnRasterSmoothingOn  = True
  res@cnFillPalette        = "BlueRed"     ; set color map
  res@lbLabelBarOn         = True        ; turn on individual lb's
  

  res@gsnStringFontHeightF = 0.015
  res@gsnLeftString  = "Leading EOF (2001-2022)"
  res@gsnRightString = sprintf("Var = "+"%5.1f", eof@pcvar) +"%"
  plot =gsn_csm_contour_map(wks,eof(0,:,:),res)
	 
  draw(plot)
  frame(wks)