theta                  = a->POT_GDS5_HYBL_10
  z                     = a->HGT_GDS5_HYBL_10
  p                     = a->PRES_GDS5_HYBL_10
  levels                = a->lv_HYBL3
  lat                   = a->g5_lat_1
  lon                   = a->g5_lon_2

 pls                    = (/ 1000.,925.,850.,800.,750.,700.,650.,600., \
                                500.,450.,400.,350.,300.,250.,200.,150.,100.,70.,50. /)

  print                 ("Calculating gradient")

T_grad_lon = theta
T_grad_lat = theta

gradsg (theta, T_grad_lon, T_grad_lat)

T_grad_lon@long_name = "longitudinal gradient (derivative)"
T_grad_lat@long_name = "latitudinal gradient (derivative)"
T_grad_lat@units = "K/m"
T_grad_lon@units = "K/m"
igradsg (T_grad_lon, T_grad_lat, theta)


print                   ("Interpolating...")
 
 T_grad_lon_int                = wrf_user_intrp3d(T_grad_lon,p,"h",pls,0.,False)
  T_grad_lat_int               = wrf_user_intrp3d(T_grad_lat,p,"h",pls,0.,False)


print ("calculate the advection")

i = 10

  u                     = f->U_GRD_GDS5_ISBL_10(i,:,:)
  v                     = f->V_GRD_GDS5_ISBL_10(i,:,:)


  lvls                  = f->lv_ISBL2

 t_lat = a->T_grad_lat_int(i,:,:)
 t_lon = a->T_grad_lon_int(i,:,:)

var = u
var = -999
var = -(v*t_lat+u*t_lon)*100
var@long_name = "heat transfer"
printVarSummary (var)