;*****************************************************
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"
;*****************************************************
; **** modify here *****
xMinF   = 694
xMaxF   = 695
varname = "amoc265"
system("ls data/micom_hm_???_??.nc >>filelist")

begin
  if (isfilepresent("amoc.nc")) then
      fid    = addfile("amoc.nc","r")
      time   = fid->time
      amoc   = fid->$varname$
  else
      files     = systemfunc("cat filelist")
      len       = dimsizes(files)
      amoc2050  = new(len,"double")
      amoc2060  = new(len,"double")
      amoc265   = new(len,"double")
      amoc450   = new(len,"double")
      amoc30S   = new(len,"double")
      time      = new(len,"double")
      amoc2050!0= "time"
      amoc2060!0= "time"
      amoc265!0 = "time"
      amoc450!0 = "time"
      amoc30S!0 = "time"
      time!0    = "time"
      fid       = addfile(files(0),"r")
      lat       = fid->lat
      latid20   = minind(abs(lat-20))
      latid50   = minind(abs(lat-50))
      latid60   = minind(abs(lat-60))
      latid265  = minind(abs(lat-26.5))
      latid450  = minind(abs(lat-45.0))
      latid30S  = minind(abs(lat-(-30.0)))
          do k=0,len-1
              print(files(k))
              fid  = addfile(files(k),"r")  
              time(k)   = fid->time
              mmflxl = fid->mmflxl(0,0,:,:)
              amoc2050(k)   = max(max(mmflxl(:,latid20:latid50)))
              amoc2060(k)   = max(max(mmflxl(:,latid20:latid60)))
              amoc265(k)    = max(max(mmflxl(:,latid265)))
              amoc450(k)    = max(max(mmflxl(:,latid450)))
              amoc30S(k)    = max(max(mmflxl(:,latid30S)))
          end do
      amoc2050=amoc2050/1e9    ; kg/s->Sv/s
      amoc2060=amoc2060/1e9    ; kg/s->Sv/s
      amoc265=amoc265/1e9    ; kg/s->Sv/s
      amoc450=amoc450/1e9    ; kg/s->Sv/s
      amoc30S=amoc30S/1e9    ; kg/s->Sv/s
      amoc2050@long_name    = "Maximum  Atlantic Meridional Overturning Circulation between 20N and 50N"
      amoc2060@long_name    = "Maximum  Atlantic Meridional Overturning Circulation between 20N and 60N"
      amoc265@long_name     = "Maximum  Atlantic Meridional Overturning Circulation at 26.5N"
      amoc450@long_name     = "Maximum  Atlantic Meridional Overturning Circulation at 45.0N"
      amoc30S@long_name     = "Maximum  Atlantic Meridional Overturning Circulation at 30.0S"
      amoc2050@units         = "Sv"
      amoc2060@units         = amoc2050@units
      amoc265@units          = amoc2050@units
      amoc450@units          = amoc2050@units
      amoc30S@units          = amoc2050@units

      time&time = time
      amoc2050&time  = time
      amoc2060&time  = time
      amoc265&time   = time
      amoc450&time   = time
      amoc30S&time   = time

      fid   = addfile("amoc.nc","c")
      fid->amoc2050 = amoc2050
      fid->amoc2060 = amoc2060
      fid->amoc265  = amoc265
      fid->amoc450  = amoc450
      fid->amoc30S  = amoc30S
      fid->time     = time

      fid   = addfile("amoc.nc","r")
      amoc  = fid->$varname$
      printVarSummary(amoc)
  end if
      yrs   = cd_calendar(time,4);
      yrs@units = "Year"

      wks_type = "pdf"
;      wks_type = "png"
;      wks_type@wkWidth = 1024
;      wks_type@wkHeight = 512
      wks  = gsn_open_wks(wks_type,"amoc")

     res                    = True

     res@gsnPaperOrientation= "Landscape"
     res@vpWidthF           = 0.6
     res@vpHeightF          = 0.6*0.618     ; Changes the aspect ratio

;     res@tiMainString               = "Column-integrated tracer content (log10)"
;     res@tiMainFontThicknessF       = 2.0 
     res@gsnStringFont              = "helvetica-bold"
     res@gsnStringFontHeightF       = 0.005    ;0.025
     res@gsnLeftString              = amoc@long_name
     res@gsnRightString             = amoc@units
;     res@gsnLeftStringFontHeightF   = 0.02
;     res@gsnRightStringFontHeightF  = 0.02
     res@tiXAxisString              = yrs@units
     res@tiYAxisString              = amoc@units

     if (any(ismissing((/xMinF,xMaxF/))) ) then
         res@trXMinF    = tointeger(min(yrs))
         res@trXMaxF    = tointeger(max(yrs))
     else
         res@trXMinF    = xMinF
         res@trXMaxF    = xMaxF
     end if

     plot    = gsn_csm_xy(wks,yrs,amoc,res)
end