;this script constructs the sst forcing data based on c4_hist2 simulation.
;the observed SST anomalies are used over the equtorial eastern ATLific Ocean;
;the buffer zone is set and the simulated sst is used over other regions 
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 "$NCARG_ROOT/lib/ncarg/nclscripts/csm/popRemap.ncl"
;---------------------------------------------------------
 begin

   system("/bin/cp /glade/work/leishanj/CESM1_2_2/sst_flow/*.nc ./")
   lat5 =  -35.5;6
   lat6 =   35.5 ;30
   lon5 =  170.5 ;30.5 ;289
   lon6 =  290.5 ;385.5 ;
   
;construct a restoring coefficient
   X2o=new((/180,360/),"float")
   X2o=0.
   X2o!0="lat"
   X2o!1="lon"
   X2o&lat=fspan(-89.5,89.5,180)
;   X2o&lon=fspan(0.5,359.5,360)
   X2o&lon=fspan(30.5,389.5,360)
   X2o&lat@units="degrees_north"
   X2o&lat@long_name="latitude"
   X2o&lon@units="degrees_east"
   X2o&lon@long_name="longitude"
   
   lat = X2o&lat
   lon = X2o&lon

	  
   X2o({lat|:},{lon|:})=0.

;the grid number for the buffer zone
   nlon  = 5
   nlat  = 5

   do ilon = 1 , nlon+1
      print(1.-ilon*1./(nlon+1))
      X2o({lat|lat5+ilon-nlon-1:lat6-ilon+nlon+1},{lon|lon5+ilon-nlon-1:lon6-ilon+nlon+1})=ilon*1./(nlon+1)
   end do
  

   X2o({lat|16:24},{lon|265:300})=0.0
   X2o({lat|9:16},{lon|275:300})=0.0
   
   X2o = X2o*10.
   printVarSummary(X2o)
   X2o = lonPivot(X2o,360.5)
   printVarSummary(X2o)
     
   X2=(/PopLatLon(X2o,"1x1d","gx1v6","bilin","da","100716")/)
   printVarSummary(X2)
	
   ff=addfile("./ru_ep.nc", "r")
   rtau=ff->rtau
   copy_VarCoords(rtau, X2)
  
   system("/bin/rm -f ./PAC_forcing.nc")
   nf=addfile("./PAC_forcing.nc", "c")
   nf->X2=X2

; -------------------------- plot -------------------------------
   wks=gsn_open_wks("pdf","./PAC_ideal_nudging-weight")
   ;gsn_define_colormap(wks,"rainbow")
   gsn_define_colormap(wks,"rainbow+white")
   
   res1 = True
   res1@gsnDraw               = False
   res1@gsnFrame              = False
   res1@gsnAddCyclic        = True
   ;res1@gsnMaximize=True
   res1@gsnPaperOrientation = "portrait"
   res1@gsnContourNegLineDashPattern=2
   res1@gsnRightString      = ""
   res1@gsnLeftStringFontHeightF=0.015
   res1@gsnRightStringFontHeightF=0.015
   res1@tiMainString        = ""
   res1@mpFillOn            = False
   res1@mpCenterLonF        = 280.
   res1@tmXBLabelFontHeightF= 0.015
   res1@tmYLLabelFontHeightF= 0.015
   ;res1@tmYLLabelFontThicknessF=0.8
   ;res1@tmXBLabelFontThicknessF=0.8
   res1@pmLabelBarOrthogonalPosF = 0.2
   ; res1@tmYLLabelStride     = 0
   res1@tmXBLabelFont       = 25
   res1@tmYLLabelFont       = 25
   ;res1@tmXBLabelDeltaF=-0.8
   ;res1@tmYLLabelDeltaF=-0.8
   ; res1@tmXTOn              = False
   ; res1@tmYROn              = False
   ; res1@tmXBAutoPrecision   = False
   ; res1@tmXBMinorOn         = False
   ; res1@tmYLMinorOn         = False
   res1@cnFillOn            = True
   res1@cnLinesOn           = False
   res1@cnLineLabelsOn      = False
   res1@cnInfoLabelOn       = False
   res1@cnLineLabelFormat   ="@^sg"
   
   res1@cnLevelSelectionMode="ExplicitLevels"
   res1@cnLevels            = (/0.2,0.6,1.0,1.4,1.8/) ;(/0.1,0.2,0.4,0.6,0.8/) ;
   res1@cnFillColors        = (/238,50,82,130,178,230/)
   res1@gsnLeftString       = "Weight"
   res1@gsnRightString      = ""
   
   lres                     = True
   lres@gsLineColor         = "black"
   lres@gsLineThicknessF    = 2.0
   lres@gsLineDashPattern   =  2
   
   plot = gsn_csm_contour_map_ce(wks,X2o,res1)
   line = gsn_add_polyline(wks,plot,(/0,180,360/),(/0,0,0/),lres) 
   draw(plot)
end