;-----------------------------------------------------------------
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" ;BRH

begin
;-- read the data and define

   spec_units = ("O3 ppb")
   factor = (/1.e+09/)
   wkdir  = "./"
   diri   = "./"
   fili   = systemfunc("ls "+diri+"f.e20.FSCD_FullEmissNig.cam.h0.2013*.nc")
   fili1  = systemfunc("ls "+diri+"f.e20.FSCD_ZeroEmissOverNig.cam.h0.2013*.nc")
   print(fili)
   print(fili1)
  minLat   = -35
  maxLat   =  35
  minLon   = -30
  maxLon   =  60

  f        = addfiles(fili, "r")
  f1       = addfiles(fili1,"r")

; calling variable

     o3x_surf = f[:]->O3(:,55,:,:)     ; input only surface
         o3x_surf1=f1[:]->O3(:,55,:,:)
   o3x_surf = o3x_surf*1.e+09        ; units change
   o3x_surf1 = o3x_surf1*1.e+09
;;o3x_surf@units = "..."
printVarSummary(o3x_surf)
printVarSummary(o3x_surf1)

; calling surface
; interpolate ozone on pressure levels
   lev1 = f[:]->lev   ;print, lev and print lev
   lat1 = f[:]->lat
   lon1 = f[:]->lon
   hyam=f[:]->hyam
   hybm=f[:]->hybm
   ps1=f[:]->PS
   p0 = 1.e+05
 levn = (/"Surface","525hpa","200hpa"/)
 ln = (/55,33,25/)

; o3x_surf = f[:]->O3(0,55,:,:)
; O3= vinth2p(O3,hyam,hybm,lev1,ps1,1,p0,1,False)

; o3x_surf1 = f[:]->O3(0,55,:,:)
; O3= vinth2p(O3,hyam,hybm,lev1,ps1,1,p0,1,False)


; synchronising dateread each file date

   date = f[:]->date
   date1 = f1[:]->date

   print(date)
   print(date1)
   printVarSummary(date)
   printVarSummary(date1)

; reorder longitudes
  o3x_surf := lonFlip(o3x_surf)
  o3x_surf1 := lonFlip(o3x_surf)

  printVarSummary(o3x_surf)
  printVarSummary(o3x_surf1)

; Simple/Crude climatological seasonal mean from a 12-month climatology
  xSea_surf = o3x_surf(0:11:3,:,:)
  xSea_surf1 = o3x_surf1(0:11:3,:,:)
  printVarSummary(xSea_surf)
  printVarSummary(xSea_surf1)

  xSea_surf(0,:,:) = (o3x_surf(0,:,:) + o3x_surf(1,:,:) + o3x_surf(11,:,:) )/3.0
  xSea_surf(1,:,:) = (o3x_surf(2,:,:) + o3x_surf(3,:,:) + o3x_surf(4,:,:) )/3.0
  xSea_surf(2,:,:) = (o3x_surf(5,:,:) + o3x_surf(6,:,:) + o3x_surf(7,:,:) )/3.0
  xSea_surf(3,:,:) = (o3x_surf(8,:,:) + o3x_surf(9,:,:) + o3x_surf(10,:,:) )/3.0

  xSea_surf1(0,:,:) = (o3x_surf1(0,:,:) + o3x_surf1(1,:,:) + o3x_surf1(11,:,:) )/3.0
  xSea_surf1(1,:,:) = (o3x_surf1(2,:,:) + o3x_surf1(3,:,:) + o3x_surf1(4,:,:) )/3.0
  xSea_surf1(2,:,:) = (o3x_surf1(5,:,:) + o3x_surf1(6,:,:) + o3x_surf1(7,:,:) )/3.0
  xSea_surf1(3,:,:) = (o3x_surf1(8,:,:) + o3x_surf1(9,:,:) + o3x_surf1(10,:,:) )/3.0


  printMinMax(xSea_surf(0,:,:), True)
  printMinMax(xSea_surf1(0,:,:), True)

; this just picks January 2013
; for annual mean you can just average over January to December

;  o3aa_surf = dim_avg_n_Wrap(O3(date_dice_i:date_dice_i+10,:,:),0)
 ; o3a_surf = dim_avg_n_Wrap(O3(date_cmip_i:date_cmip_i+10,:,:),0)

  O3diff0 = xSea_surf(0,:,:)
  o3diff0 = xSea_surf(0,:,:)-xSea_surf1(0,:,:)

  O3diff1 = xSea_surf(1,:,:)
  o3diff1 = xSea_surf(1,:,:)-xSea_surf1(1,:,:)

  O3diff2 = xSea_surf(2,:,:)
  o3diff2 = xSea_surf(2,:,:)-xSea_surf1(2,:,:)

  O3diff3 = xSea_surf(3,:,:)
  o3diff3 = xSea_surf(3,:,:)-xSea_surf1(3,:,:)

  printMinMax(O3diff0,0)
 printVarSummary(O3diff0)

 printMinMax(O3diff1,0)
 printVarSummary(O3diff1)

        factor=1.e+10
        o3x_surf = o3x_surf*factor

        factor=1.e+10
        o3x_surf1 = o3x_surf1*factor

;  the workstation (plot type and name)
  wks = gsn_open_wks("x11", "O3_FullEmissNig_vs_ZeroEmissOverNig")

;-- set resources

    res                 = True
 ;gsn_merge_colormaps(wks,"amwg","amwg_blueyellowred") ; merge two color maps
 gsn_define_colormap(wks,"BlueWhiteOrangeRed")
 res@gsnDraw                = False
 res@gsnFrame               = False
 res@gsnMaximize = False
  res@cnFillOn             = True                  ; turn on color fill
  res@cnLinesOn         = False                 ; no contour lines
  res@cnLineLabelsOn    = False
  res@tmXTOn    = False
  res@tmXBOn = False
  res@tmYLLabelFontHeightF = 0.015
  res@tmYRMode             = "Automatic"          ; turn off special labels on right axis
  res@lbLabelBarOn      = True
  res@lbOrientation     = "Vertical"
  res@cnInfoLabelOn        = False
  res@cnLevelSelectionMode      = "ExplicitLevels"       ; manual levels
  res@cnLevelSelectionMode  = "ManualLevels"   ;-- set contour levels manually
  res@cnMinLevelValF        =-145            ;-- minimum contour level
  res@cnMaxLevelValF        =  95
  res@cnMissingValFillColor = 0
  res@txFontHeightF = 0.020
  res@txFont = "helvetica-bold"
  res@gsnPaperOrientation = "landscape"
  res@pmLegendWidthF = 0.15
  res@pmLegendHeightF = 0.15
  res@lgLabelFontHeightF = .015
  res@lgPerimOn = True
  res@txFontHeightF = 0.015
  res@mpDataSetName         = "Earth..4"   ; This new database contains
                                           ; divisions for other countries.
  res@mpDataBaseVersion     = "MediumRes"  ; Medium resolution database

  res@mpOutlineOn           = True         ; Turn on map outlines
  res@mpOutlineSpecifiers   = (/"Africa:states"/)
  res@mpLimitMode           = "LatLon"
  res@mpMaxLatF             =  35          ;  Africa limits
  res@mpMinLatF             = -35
  res@mpMaxLonF             = 60
  res@mpMinLonF             = -30

;*********************************************************************
   res@cnLevelSelectionMode = "ExplicitLevels"
   res@cnFillOn             = True
   res@cnLinesOn            = False
   res@cnLineLabelsOn       = False
   res@lbTitleOn = True
   res@lbTitleString = "ppb"
   res@lbLabelFontHeightF = 0.018
   res@lbTitleFontHeightF = 0.02
   res@lbLabelStride  = 2
   res@lbOrientation = "Horizontal"
; contour levels
 con = new(14,float)
 con=(/0.5, 1.,2., 3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13./)
 diffcon = new(14,float)
; diffcon=(/-140.,-130., -100.,-80.,-60.,-20,5,20.,30.,40., 50.,60.,70.,75./)
  diffcon=5*(/-6.,-5., -4.,-3.,-2.,-1.,-0.5,0.5,1.,2., 3.,4.,5.,6./)

  res@cnLevels = con
  res@mpOutlineOn           = True         ; Turn on map outlines
  res@mpOutlineSpecifiers   = (/"Africa:states"/)
  res@mpLimitMode           = "LatLon"
  res@mpMaxLatF             =  35          ;  Africa limits
  res@mpMinLatF             = -35
  res@mpMaxLonF             = 60
  res@mpMinLonF             = -30

;**********************************************************************

  res@cnLevels = diffcon

 plot0 = gsn_csm_contour_map_ce(wks,O3diff0,res)
         res@gsnRightString = "O3 ppb"
         res@gsnLeftString= "DJF_dry cold"
         res@gsnCenterString = "Full minus Zero Nigeria"
         res@txFontHeightF   = 0.011
 plot1 = gsn_csm_contour_map_ce(wks,O3diff1,res)
          res@gsnRightString = "O3 ppb"
          res@gsnLeftString= "MAM_dry warm"
          res@gsnCenterString = "Full minus Zero Nigeria"
          res@txFontHeightF   = 0.011
 plot2 = gsn_csm_contour_map_ce(wks,O3diff1,res)
          res@gsnRightString = "O3 ppb"
          res@gsnLeftString= "JJA_monsoon"
          res@gsnCenterString = "Full minus Zero Nigeria"
          res@txFontHeightF   = 0.011
 plot3 = gsn_csm_contour_map_ce(wks,O3diff1,res)
          res@gsnRightString = "O3 ppb"
          res@gsnLeftString= "SON_post_monsoon"
          res@gsnCenterString = "Full minus Zero Nigeria"
          res@txFontHeightF   = 0.011
;--Panel all three plots in a m x n configuration
pres = True

; add lat long info to plot
                ; add lat long info to plot
    txres               = True
    txres@txFontHeightF = 0.017

   panres = True
   panres@gsnFrame = False
   panres@gsnMaximize = True
    panres@gsnPanelTop = 0.96
     panres@gsnPanelBottom = 0.05
     gsn_panel (wks, (/plot0,plot1,plot2,plot3/),(/2,2/),panres)

    frame(wks)
;-------------------------------------------------------------------------
 delete(res)


end