datapath = "./"
varnames = (/"AODVIS","colccn.3"/)
nvar     = dimsizes(varnames)

filenames = (/"v2_ndg_cdnc_ssat_diag_simple_pi.eam.h0.2011-01.nc", \
              "v2_ndg_cdnc_ssat_diag_simple_pd.eam.h0.2011-01.nc"/)
nfile = dimsizes(filenames)
Files = new( nfile,"file" )

nPlotColumn = 5
plot = new( (/nvar,nPlotColumn/),"graphic")

wks = gsn_open_wks("pdf","./compare_two_runs")

res = True
res@gsnDraw  = False
res@gsnFrame = False
res@cnFillOn = True
res@cnLinesOn = False
res@mpProjection           = "WinkelTripel"
res@mpOutlineOn            = True
res@mpPerimOn              = False
res@mpGridAndLimbOn        = True
res@mpGridLineDashPattern  = 2
res@gsnStringFontHeightF   = 0.02

;----------------------------------------------------------
; Get lat/lon info from the first file. Also check if the
; two files have the same horizontal dimension size.
;----------------------------------------------------------
do ifile = 0,nfile-1

   File = addfile(datapath+"/"+filenames(ifile),"r")

   if (ifile.eq.0) then

      area = File->area          ; read in the "area" var
      ncol = dimsizes(area)      ; get the number of grid columns
      res@sfXArray  = File->lon  ; get longitude info needed for plotting
      res@sfYArray  = File->lat  ; get latitude  info needed for plotting

   else

      ncol1 = dimsizes(File->area)  ; the the number of grid columns

      if (ncol1.ne.ncol) 
         print("Error: two files have different numbers of grid columns. Abort.")
         exit
      end if

   end if

end do

;----------------------------------------------------------
; Read variables. Calculate and show differences.
;----------------------------------------------------------
FillValue = -999. ; for avoiding division by zero

do ivar = 0,nvar-1

   array4plotting = new( (/nPlotColumn,ncol/),"float" )
   res@gsnLeftString = "" 
   res@gsnRightString = ""
   res@tiMainString = varnames(ivar) 


   res@cnFillPalette = "WhViBlGrYeOrRe"
   do ifile = 0,nfile-1
      File = addfile(datapath+"/"+filenames(ifile),"r")
      array4plotting(ifile,:) = File->$varnames(ivar)$

      if (ifile.eq.0) then
         res@gsnCenterString = "Ctrl"
      else
         res@gsnCenterString = "Test"
      end if
      plot(ivar,ifile) = gsn_csm_contour_map(wks,array4plotting(ifile,:),res)
   end do

   res@cnFillPalette = "BlueDarkRed18"

   ; difference

   ic = 2

   res@gsnCenterString = "Difference: test - ctrl"
   array4plotting(ic,:) = array4plotting(1,:) - array4plotting(0,:)
   symMinMaxPlt(array4plotting(ic,:),14,False,res)

   plot(ivar,ic) = gsn_csm_contour_map(wks,array4plotting(ic,:),res)

   ; relative difference (1)

   ic = 3
   res@gsnCenterString = "Relative difference: (test - ctrl)/ctrl"

   tmp = array4plotting(0,:)                  ; use the first simulation for normalization
   denom = where( tmp.ne.0., tmp, FillValue ) ; mask out zeros
   denom@_FillValue = FillValue

   array4plotting(ic,:) = ( array4plotting(1,:) - array4plotting(0,:) )/denom

   symMinMaxPlt(array4plotting(ic,:),14,False,res)
   plot(ivar,ic) = gsn_csm_contour_map(wks,array4plotting(ic,:),res)

   ; relative difference (2)

   ic = 4
   res@gsnCenterString = "Relative difference: (test - ctrl)*2/(test + ctrl)"

   tmp = 0.5*( array4plotting(0,:) + array4plotting(1,:) )  ; use the average of the two runs for normalization
   denom = where( tmp.ne.0., tmp, FillValue )               ; mask out zeros
   denom@_FillValue = FillValue

   array4plotting(ic,:) = ( array4plotting(1,:) - array4plotting(0,:) )/denom

   symMinMaxPlt(array4plotting(ic,:),14,False,res)
   plot(ivar,ic) = gsn_csm_contour_map(wks,array4plotting(ic,:),res)

   ; clean up

   delete(res@cnLevelSelectionMode)
   delete(res@cnMinLevelValF)
   delete(res@cnMaxLevelValF)
   delete(res@cnLevelSpacingF)
   delete(tmp)
   delete(denom)
   delete(array4plotting)
end do

resp = True
resp@gsnMaximize = True
resp@gsnPanelXWhiteSpacePercent = 2
resp@gsnPanelYWhiteSpacePercent = 5
gsn_panel(wks,ndtooned(plot),(/nvar,nPlotColumn/),resp)