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/cnmap/cnmap.ncl"

begin
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;input
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
PDataSST=addfile("~/data/HadISST_sst.nc","r")
sst=PDataSST->sst
target_month=8
sst!1="lat"
sst!2="lon"

latmin=-30
latmax=40
lonmin=30
lonmax=135

sst_for_eof_1=sst(lat|:,lon|:,time|1307+target_month:1739:12)
sst_for_eof_1_clim=clmMonLLT(sst_for_eof_1)
sst_for_eof_2=calcMonAnomLLT(sst_for_eof_1,sst_for_eof_1_clim)
sst_for_eof_dtrend=dtrend_n(sst_for_eof_2,False,2)
copy_VarCoords(sst_for_eof_1,sst_for_eof_dtrend)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;process
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;pre-process
; =================================================================
; create weights:  sqrt(cos(lat))   [or sqrt(gw) ]
;                                     =                                                                       =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =
  rad                                 = 4.*atan(1.)/180.
  clat                                = PDataSST->latitude
  clat                                = sqrt( cos(rad*clat) )                 ; gw for gaussian grid
  clat!0                              = "lat"
sst_for_eof_dtrend_weighted           = sst_for_eof_dtrend                                   ; copy meta data
sst_for_eof_dtrend_weighted           = sst_for_eof_dtrend*conform(sst_for_eof_dtrend, clat, 0)
;sst_for_eof_dtrend_weighted@long_name = "Wgt: "+sst_for_eof_dtrend_weighted@long_name
sst_for_eof=sst_for_eof_dtrend_weighted({lat|latmin:latmax},{lon|lonmin:lonmax},time|:)


neof=3
optEOF=True
;optEOF@jopt=1
sst_eof_v=eofunc_Wrap(sst_for_eof,neof,optEOF)
sst_eof_ts=eofunc_ts_Wrap(sst_for_eof,sst_eof_v,False)
printVarSummary(sst_eof_v)
printVarSummary(sst_eof_ts)
asciiwrite("eof_ts_1_"+target_month+".txt",sst_eof_ts(0,:))

stddev_eof_ts_1=dim_stddev_Wrap(sst_eof_ts(0,:))
print(stddev_eof_ts_1)
sst_rebuild_1=sst_eof_v*stddev_eof_ts_1*2.5
copy_VarCoords(sst_eof_v,sst_rebuild_1)
sst_rebuild=sst_rebuild_1(0,:,:)
printVarSummary(sst_rebuild)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;resourse
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
res=True
res@gsnFrame=False
res@gsnDraw=False
res@gsnAddCyclic         = False        ; data not cyclic
res@gsnSpreadColors      = True
;res@gsnSpreadColorStart=-1
;res@gsnSpreadColorEnd=2
res@cnFillOn=True
res@cnLinesOn=False
res@cnLineLabelsOn=False
res@cnLevelSelectionMode="ExplicitLevels"
res@cnLevels=fspan(-0.8,0.8,17)
res@mpMinLatF             = latmin
res@mpMaxLatF             = latmax
res@mpMinLonF             = lonmin
res@mpMaxLonF             = lonmax
res@mpDataSetName         = "Earth..4"   ; This new database contains
res@mpDataBaseVersion     = "MediumRes"  ; Medium resolution database
res@mpOutlineOn           = True         ; Turn on map outlines
res@mpOutlineSpecifiers   = (/"China:states","Taiwan"/)       ;China:states
res@gsnLeftString=""
res@lbLabelBarOn=True

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;plot
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
wks=gsn_open_wks("png",get_script_prefix_name+"_"+target_month)
;gsn_define_colormap(wks,"matlab_jet")
res@gsnRightString = sprintf("%5.1f", sst_eof_v@pcvar(0)) +"%"
;plot=gsn_csm_contour_map(wks,sst_eof_v(0,:,:),res)
plot=gsn_csm_contour_map(wks,sst_rebuild,res)

draw(plot)
frame(wks)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;add to amip clm
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
PInputData=addfile("./amip2_sst.data.climo.nc","r")
;system ("rm -f "+"~/FAMIL/exp/iobm/input/ocn/sst/amip2_sst.data.climo.nc")
;POutputData=addfile("~/FAMIL/exp/iobm/input/ocn/sst/amip2_sst.data.climo.nc","c")
system ("rm -f "+"./amip2_sst.data.climo.nc")
POutputData=addfile("./amip2_sst.data.climo.nc","c")

;weights matrix
weight_dim=dimsizes(sst_rebuild)
weight=new(weight_dim,float)
weight=0
bond_wid=5
do i=1,bond_wid
weight(i:weight_dim(0)-i,i:weight_dim(1)-i)=weight(i:weight_dim(0)-i,i:weight_dim(1)-i)+1./bond_wid
end do
copy_VarCoords(sst_rebuild,weight)
sst_rebuild=sst_rebuild*weight


sst_mod=short2flt(PInputData->sst)
printVarSummary(sst_mod(time|:,{lat|latmin:latmax},{lon|lonmin:lonmax}))
sst_mod(time|target_month-1,{lat|latmin:latmax},{lon|lonmin:lonmax})=sst_mod(time|target_month-1,{lat|latmin:latmax},{lon|lonmin:lonmax})+sst_rebuild

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;output
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
dim_names = (/"time","lat","lon"/)
;time dimension is unlimited, but currently set to 1.
dim_sizes = dimsizes(sst)
dim_unlimited = (/True,False,False/)

print("Defining dimension sizes to be " + dim_sizes)
filedimdef(POutputData,dim_names,dim_sizes,dim_unlimited)


  att_names = getvaratts(PInputData)           ; get CCM file's global attributes
  if(.not.all(ismissing(att_names))) then
    do i = 0,dimsizes(att_names)-1
      print("copy_fileatts: global attributes->" + att_names(i))
      POutputData@$att_names(i)$ = PInputData@$att_names(i)$     ; copy CCM file's global att ributes
    end do
  end if

  names = getfilevarnames(PInputData)       ; get CCM file's variable names and print them
  print (names)
;
; loop over variables and copy to netCDF file
;
  do i=0, 2
    print("writing POutputData: i,name="+ i+ "   "+ names(i))
    POutputData->$names(i)$ = PInputData->$names(i)$
  end do
    POutputData->sst= sst_mod
  do i=4, dimsizes(names)-1
    print("writing POutputData: i,name="+ i+ "   "+ names(i))
    POutputData->$names(i)$ = PInputData->$names(i)$
  end do

  print (POutputData)      ; print overview of contents of netCDF file
end