;---------------------------------------------------------------------------
;---------------------------------------------------------------------------
;   Script to output ERSST for CPC nino-3.4 index         ;
;   and power spectra.			                  ;
;----------------------------------------------------------
; There are four region based indices used to monitor the tropical Pacific: 
; Nino 1+2 (0-10S, 90W-80W), Nino 3 (5N-5S, 150W-90W), 
; Nino 3.4/ONI (5N-5S, 170W-120W) and Nino 4 (5N-5S, 160E-150W)
;
; NOAA's operational definitions of El Nino and La Nina conditions are based
; upon the Oceanic Nino Index [ONI]. The ONI is defined as the 3-month running
; means of SST anomalies in the Nino 3.4 region [5N-5S, 120-170W]. The anomalies
; are derived from the 1971-2000 SST climatology.
;---------------------------------------------------------------------------
	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
	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"
;---------------------------------------------------------------------------
	begin
;
  	latS   = -62.0
  	latN   =  62.0
  	lonL   = 30.0
  	lonR   = 300.0
  
  	nrun   = 3             ; length of running average

	run1   = getenv("MODE")
        ymStrt = (stringtointeger(getenv("ISYY")))*100+1
        ymLast = (stringtointeger(getenv("IEYY")))*100+12
        isy    = ymStrt/100 
        iey    = ymLast/100

	if(run1.eq."historical") then
           dir = "/archive/CMIP6/TaiESM1/historical/atmos/mon/r1i1p1f1/"
        else
;          dir = "/lfs/archive/taiesm_ppost/"+run1+"/v20210908/"
           dir = getenv("path")
        end if
   dir = getenv("path")
   model = getenv("model")

;---------------------------------------------------------------------
        ifn    = systemfunc("ls "+dir+"ts_Amon_*.nc")
        a      = addfiles(ifn,"r")
        data   = a[:]->ts;(:,{latS:latN},{lonL:lonR})
;---------------------------------------------------------------------
  	yyyymm = ut_calendar(data&time,-1)
  	yyyymm!0 = "date"
  	yyyymm&date = yyyymm
 	delete(data&time)
  	data&time = yyyymm

  	tStrt = ind(yyyymm.eq.ymStrt)      ; indices of selected times
  	tLast = ind(yyyymm.eq.ymLast)

  	x = data(tStrt:tLast,{latS:latN},{lonL:lonR})
        x = x - 273.15
        x@units = "degC"
;*********************************
; 	Climatology and anomalies from base climatology   
;*********************************
   	xClm     = clmMonTLL(x) ;(iClmStrt:iClmLast,{latS:latN},{lonL:lonR}))

        printVarSummary(xClm)

  	xAnom    = calcMonAnomTLL (x,  xClm ) 
  	xAnom@long_name = "SST Anomalies"
  	;rintVarSummary(xAnom)

;*********************************
; 	de-trend
;*********************************

   	dsst = dtrend_leftdim(xAnom,False)
   	copy_VarCoords(xAnom,dsst)
  	;rintVarSummary(dsst)

;       dsst = xAnom
;*********************************
; 	cal. STD
;*********************************

        std = dim_stddev_n_Wrap(dsst,0)
        printVarSummary(std)

;	maks land

 	idd  = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc","r")
   	lsdata = idd->LSMASK
   	lsm  = landsea_mask(lsdata,std&lat,std&lon)
   	std  = mask( std,lsm.eq.1,False)
   	xClm = mask(xClm,lsm.eq.1,False)
;	output data

        ofn = model+"_"+run1+"_map-STD_"+ymStrt+"-"+ymLast+".nc"
        system("rm -rf "+ofn)
        odd = addfile(ofn,"c")
        odd->std = std
	odd->xClm = xClm
;*********************************
;       setting plot resoure
;*********************************

       wks     = gsn_open_wks("pdf",model+"_"+run1+"_map-STD_"+ymStrt+"-"+ymLast)

        res                = True
        res@gsnDraw        = False        ; don't draw yet
        res@gsnFrame       = False        ; don't advance frame yet

        res@gsnAddCyclic   = False

        res@cnLinesOn      = False
        res@cnLineLabelsOn = False
        res@cnFillOn       = True
        res@lbLabelsOn     = True
        res@lbLabelBarOn   = True ;False

        res@tmXTOn       = False
;       res@tmYROn      = False
        res@tmYLMinorOn   = False
        res@tmXBMinorOn   = False

        res@gsnStringFontHeightF = 0.020
        res@tmXBLabelFontHeightF = 0.015
        res@tmYLLabelFontHeightF = 0.015

        res@cnFillPalette        = "precip3_16lev"

;        res@cnLevelSelectionMode = "ExplicitLevels"   ; set explicit contour levels
;        res@cnLevels             = (/-0.5,-.4,-.3,-.2,-.1,.1,.2,.3,.4,.5/)

        res@cnLevelSelectionMode  = "ManualLevels"
        res@cnMinLevelValF        = 0.2    
        res@cnMaxLevelValF        =  1.6    
        res@cnLevelSpacingF       =  .2

        res@mpCenterLonF         = 180
        res@mpMinLatF            = latS
        res@mpMaxLatF            = latN
        res@mpMinLonF            = lonL
        res@mpMaxLonF            = lonR

        res@mpFillOn             = True 

;       set plot psi contour res

        res2                  = True
        res2@gsnDraw          = False        ; don't draw yet
        res2@gsnFrame         = False        ; don't advance frame yet

        res2@gsnAddCyclic     = False
        res2@cnLinesOn        = True
        res2@cnLineLabelsOn   = False
        res2@cnLineColor      = "white"  ;"lightskyblue4"
        res2@cnLineThicknessF = 1.8
        res2@cnInfoLabelOn    = False
        res2@cnFillOn         = False

        res2@tiMainString     = " "

        res3 = res2

        res2@cnLineLabelBackgroundColor = -1

        res2@cnLevelSelectionMode  = "ManualLevels"
        res2@cnMinLevelValF        = 3.
        res2@cnMaxLevelValF        = 33.
        res2@cnLevelSpacingF       = 3.


        res2@cnLineDashSegLenF          = 0.18          ; assist in controlling concentration

; 	set clm 27 degC contour
         
        res3@cnLineThicknessF     = 4.
        res3@cnLevelSelectionMode = "ExplicitLevels"
        res3@cnLevels     =(/27./) 
        res3@gsnLeftString = " "
        res3@gsnRightString = " "
        res3@gsnCenterString = " "

        res@gsnLeftString        = model+" "+run1+" monthly STD. ("+ymStrt/100+"-"+ymLast/100+")"
        plot = gsn_csm_contour_map_overlay(wks, std, dim_avg_n_Wrap(xClm,0),res,res2)
        cplot= gsn_csm_contour(wks,dim_avg_n_Wrap(xClm,0),res3)
      
        overlay(plot ,cplot)
        draw(plot)
        frame(wks)

        end