;TB is the total mass using function sum
;TB2 is the total mass using demensional summation first
;I'm expection Tb and TB2 to be the same values or at least close. But it is not the case...
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/wrf/WRFUserARW.ncl"


  case = "33.3m_1536"
  var = "QVAPOR"

  dir = "/glade/scratch/xchu/ideal_case/BLT/"+case+"/joint_wrf/aux_result/"
  fileTB = systemfunc("ls "+dir+"auxhist*")

  tsize = dimsizes(fileTB)

  TB = new((/tsize/),double)
  TB2 = new((/tsize/),double)

  do tt = 0,2;tsize-1
     a=addfile(fileTB(tt),"r")
     dx = a@DX
     dy = a@DY
     times  = wrf_user_list_times(a)

;get dimension
     mdims1 = getfilevardimsizes(a,"QVAPOR")
;process for dz
     phb1  = wrf_user_getvar(a, "PHB", 0)
     ph1   = wrf_user_getvar(a, "PH", 0)
     hgt_m1=(phb1+ph1)/9.81
     delete(phb1)
     delete(ph1)
     dzTB=new((/mdims1(1),mdims1(2),mdims1(3)/),"double")
     do k=0,mdims1(1)-1
      dzTB(k,:,:)=(hgt_m1(k+1,:,:)-hgt_m1(k,:,:))
     end do
     delete(hgt_m1)

     varTB=wrf_user_getvar(a,var,0)/wrf_user_getvar(a,"ALT",0)*dzTB*dx*dy 

     tempTB  = dim_sum_n(varTB,(/1,2/))
     TB(tt)  = sum(tempTB)
     TB2(tt) = sum(varTB)
  end do

  print(TB(0:2))
  print(TB2(0:2))