; Atlantic 70W-20E (103:127 and 0-6)
; Pacific 120E-70W (42-103)
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/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
  in = addfile("/ptmp/notaro/pimodTC1/ho.foam1-5lpj.pimodTC1.anni10351134.nc","r")
  v  = in->V   
  lat= in->lat
  lon= in->lon
  lev= v&lev
  time = in->time

  v!0="time"
  v!1="lev"
  v!2="lat"
  v!3="lon"
  v&time=time
  v&lev=lev
  v&lat=lat
  v&lon=lon

  nlev = dimsizes(lev)
  nlat = dimsizes(lat)
  dx = 0.313084e6
  v@_FillValue = v@missing_value  

; calculate dz array since not on file
  dz = new(nlev,typeof(lev))
  do k = 0,nlev-2
    dz(k) = lev(k+1)-lev(k)
  end do
  dz(nlev-1) = 0 
;**************************************************************
  d2rad = 0.017453         ; degrees to radians

  early=dim_avg(v(lev|:,lat|:,lon|:,time|0:9))
  late=dim_avg(v(lev|:,lat|:,lon|:,time|90:99))

;**************************************************************  
; calculate first intergral 
; int[lon1:lon2]v*cos(lat)*dx*dz
; this calculation is done on the z_t grid
;**************************************************************
  zone_int1 = new((/nlat,nlev/),typeof(v))     ; allocate space
  do k = 0, nlev-1
    do j = 0, nlat-1
      zone_int1(j,k) = dim_sum(early(k,j,42:103)*cos(lat(j)*d2rad)*dx*dz(k)) 
    end do
  end do
;************************************************************** 
; calculate second integral (partial summation) over levels on z_w grid
; psi(k,y)=int[k:0]zone_int 
;**************************************************************
  moc1      = new((/nlev,nlat/),typeof(v))      ; allocate space
  moc1(0,:) = 0.                                ; bottom is zero 
  do k=1,nlev-1
    moc1(k,:) = -1.0 * dim_sum(zone_int1(:,0:k))
  end do

;************************************************************** 
  moc1=moc1/1e6 ;     Sv=1e6 m3/s
  moc1!0         = "depth"
  moc1!1         = "lat"
  moc1&depth     = v&lev
  moc1&lat       = v&lat
  moc1@long_name = "eulerian merdional overturning"
  moc1@units     = "Sv"
  moc1@_FillValue=-999.

  print(max(moc1))
  print(min(moc1))
;*********************************

  zone_int2 = new((/nlat,nlev/),typeof(v))     ; allocate space
  do k = 0, nlev-1
    do j = 0, nlat-1
      zone_int2(j,k) = dim_sum(late(k,j,42:103)*cos(lat(j)*d2rad)*dx*dz(k))
    end do
  end do
  moc2      = new((/nlev,nlat/),typeof(v))      ; allocate space
  moc2(0,:) = 0.                                ; bottom is zero
  do k=1,nlev-1
    moc2(k,:) = -1.0 * dim_sum(zone_int2(:,0:k))
  end do
  moc2=moc2/1e6 ;     Sv=1e6 m3/s
  copy_VarCoords(moc1,moc2)
  moc2@_FillValue=-999.

  print(max(moc2))
  print(min(moc2))

;;;;;;;;;;;;;;;

  diff=new((/nlev,nlat/),typeof(v))
  do k = 0, nlev-1
    do j = 0, nlat-1
      if ( .not.ismissing(moc1(k,j))) then
        if (fabs(moc1(k,j)).gt.0.) then
          diff(k,j)=(moc2(k,j)-moc1(k,j))/fabs(moc1(k,j))*100.
        else
          diff(k,j)=-999.
        end if
      else
        diff(k,j)=-999.
      end if 
    end do
  end do
  copy_VarCoords(moc1,diff)
  diff@_FillValue=-999.
  print(max(diff))
  print(min(diff))
  diff=diff>-100.
  diff=diff<100.
  diff=smth9(diff,.5,.25,True)

  print(moc1(7,36))
  print(moc2(7,36))
  print(diff(7,36))
;;;;;;;;;;;;;
  wks = gsn_open_wks("ps","merid_overturn_pimodTC1_pac_diff")                 ; open a ps file
  plot=new(3,graphic)
  gsn_define_colormap(wks,"ViBlGrWhYeOrRe")      ; choose colormap
  res                         = True             ; plot mods desired
  res@gsnDraw = False
  res@gsnFrame = False 
  res@cnFillOn                = True             ; turn on color fill
  res@cnLineLabelsOn          = False            ; turns off contour line labels
  res@cnLinesOn               = False            ; turn off contour lines
  res@cnInfoLabelOn           = False            ; turns off contour info label
  res@gsnSpreadColors         = True             ; use full colormap
  res@lbOrientation           = "vertical"       ; vertical label bar
  res@lbLabelAutoStride       = True             ; nice label bar labels
  res@txFontHeightF      = 0.015   
  
  res@sfXArray                = moc1&lat          ; uses lat_t as plot x-axis
  res@sfYArray                = moc1&depth/1000.   ; convert cm to km
  res@tiXAxisString           = "Latitude"
  res@tiYAxisString           = "Depth(km)"
  res@cnMissingValPerimOn     = True             ; turn on perimeter
  res@cnMissingValFillPattern = 3                ; choose a fill pattern
  res@cnMissingValFillColor   = "black"          ; choose a pattern color 

  res@cnLevelSelectionMode    = "ExplicitLevels"   ; manually set contour levels
  res@cnLevels=(/-40.,-30.,-20.,-10.,0.,10.,20.,30.,40./)
  res@tiMainString = "Pacific Merid Overturning (Sv) PIMODTC1 Yr1-10"
  res@gsnCenterString=""
  res@gsnLeftString=""
  res@gsnRightString="" 
  res@trXMinF = -50.
  res@trXMaxF = 70.
  plot(0)  = gsn_csm_contour(wks,moc1,res)           ; create plot

  res@tiMainString = "Pacific Merid Overturning (Sv) PIMODTC1 Yr91-100"
  plot(1)  = gsn_csm_contour(wks,moc2,res)           ; create plot

  res@tiMainString = "% Change"
  res@cnLevels=(/-40.,-20.,-10.,-5.,0.,5.,10.,20.,40./)
  delete(res@cnMissingValFillPattern)
  plot(2)  = gsn_csm_contour(wks,diff,res)           ; create plot

  pres             = True                 ; mod panel plot
  pres@gsnDraw=True
  pres@gsnFrame=True
  pres@gsnMaximize = True                        ; blow up plot
  pres@gsnPaperWidth=8.5
  pres@gsnPaperHeight=11.0
  pres@gsnPaperOrientation="portrait"
  pres@gsnPaperMargin=0.07
  pres@gsnPanelLeft = 0.04
  pres@gsnPanelRight = 0.96
  pres@gsnPanelTop = 0.96
  pres@gsnPanelBottom = 0.04
  pres@gsnPanelYWhiteSpacePercent = 3
  pres@gsnPanelRowSpec=True
  gsn_panel(wks,plot,(/2,1/),pres)        ; create panel plot

end