[ncl-talk] Calculating zonal, meridional, and vertical gradients for WRF Output

Tue Dec 18 14:52:45 MST 2018

Hi,

I am trying to calculate zonal, meridional, and vertical gradients for WRF
output, and have encountered a problem along the way. For some reason , for
one of the variables I am calculating, the result for all dimensions is the
missing value of 9.96921e+36!

Below is my code for NCL.6.5. In order to calculate zonal and meridional
gradients, I'm using grad_latlon_cfd, and in order to be able to use this
function, the curvilinear wrf data is converted to rectilinear using
rcm2rgrid_Wrap.

 begin

  a = addfile("wrfout_d02_2017-07-17_00:00:00","r")

  type = "pdf"

  wks = gsn_open_wks(type,"Kinematics")
  times  = wrf_user_list_times(a)
  time   = 40

  nplots     = 4
  plot       = new(nplots,"graphic")
   res                      = True
   res at tiMainOn             = False
;--------------------------- Shear
-----------------------------------------;
     P   = wrf_user_getvar(a,"P",time)
     PB  = wrf_user_getvar(a,"PB",time)
     T   = wrf_user_getvar(a,"T",time)
     U   = wrf_user_getvar(a,"U",time)
     V   = wrf_user_getvar(a,"V",time)
     lon = wrf_user_getvar(a,"XLONG",time)
     lat = wrf_user_getvar(a,"XLAT",time)

     T = T + 300.
     P = P + PB

     TH = pot_temp(P, T, -1, False)
     Theta = TH(:,:,389:391)
     US = U(:,:,389:391)
     VS = V(:,1:,389:391)  ; The domain is 900*900 grid points.
     LAT = lat(:,390)
     LON = lon(574,389:391)

;------------- Convert Curvilinear to Rectilinear -----------
     Theta_rec =
rcm2rgrid_Wrap(lat(:,389:391),lon(:,389:391),Theta,LAT,LON,1)
     US_rec    = rcm2rgrid_Wrap(lat(:,389:391),lon(:,389:391),US,LAT,LON,1)
     VS_rec    = rcm2rgrid_Wrap(lat(:,389:391),lon(:,389:391),VS,LAT,LON,1)
;-------------------- Gradients -----------------------------
     Thgrad = grad_latlon_cfd(Theta_rec,LAT,LON,False,False)
     Ugrad  = grad_latlon_cfd(US_rec,LAT,LON,False,False)
     Vgrad  = grad_latlon_cfd(VS_rec,LAT,LON,False,False)
     dThdx = Thgrad[1]
     dThdy = Thgrad[0]

     deltaTh = ((dThdx ^2) + (dThdy ^2))^(-0.5)
     dUdy  = Ugrad[0]
     dVdx  = Vgrad[1]

     Sh = - deltaTh *(dThdx * dThdy*(dUdy + dVdx))
    ; The above shear term calculated for a S-N cross section centered on i
= 390
     Sh_avg = dim_avg_n(Sh,2)
;------------------Confluence ---------------------------------
     dVdy  = Vgrad[0]
     dUdx  = Ugrad[1]
     Con = - deltaTh *(((dThdy ^2) * dVdy) + ((dThdx ^2) * dUdx))
     Con_avg = dim_avg_n(Con,2)
;------------------- Tilting ------------------------
     W   = wrf_user_getvar(a,"W",time)
     Z   = wrf_user_getvar(a,"z",time)
     lev = Z(:,:,389:391)
     WS = W(:,:,389:391)
     WS_rec = rcm2rgrid_Wrap(lat(:,389:391),lon(:,389:391),WS,LAT,LON,1)
     lev_rec = rcm2rgrid_Wrap(lat(:,389:391),lon(:,389:391),lev,LAT,LON,1)

     dThdz = center_finite_diff_n(Theta_rec,lev_rec,False,0,0)

     Wgrad  = grad_latlon_cfd(WS_rec,LAT,LON,False,False)
     dWdy  = Wgrad[0]
     dWdx  = Wgrad[1]

     Tilt = - deltaTh *(dThdz * (dWdy(1:,:,:) * dThdy + dWdx(1:,:,:) *
dThdx))
     Tilt_avg = dim_avg_n(Tilt,2)
;------------------------ Diabatic Heating Term ----------------------

     P2  = wrf_user_getvar(a,"P",time-1)
     P3  = wrf_user_getvar(a,"P",time+1)

     PB2  = wrf_user_getvar(a,"PB",time-1)
     PB3  = wrf_user_getvar(a,"PB",time+1)

     lat2 = wrf_user_getvar(a,"XLAT",time-1)
     lat3 = wrf_user_getvar(a,"XLAT",time+1)

     lon2 = wrf_user_getvar(a,"XLONG",time-1)
     lon3 = wrf_user_getvar(a,"XLONG",time+1)

     T2  = wrf_user_getvar(a,"T",time-1)
     T3  = wrf_user_getvar(a,"T",time+1)

     Tr2 = T2(:,:,389:391) + 300.
     Tr3 = T3(:,:,389:391) + 300.

     Pr2 = P2(:,:,389:391) + PB2(:,:,389:391)
     Pr3 = P3(:,:,389:391) + PB3(:,:,389:391)

     levt = Z(:,:,389:391)

     ThetaT2 = pot_temp(Pr2, Tr2, -1, False)
     ThetaT3 = pot_temp(Pr3, Tr3, -1, False)
     LAT2 = lat2(:,390)
     LON2 = lon2(574,389:391)

     LAT3 = lat3(:,390)
     LON3 = lon3(574,389:391)

rcm2rgrid_Wrap(lat(:,388:392),lon(:,388:392),ThetaT1,lat(:,390),lon(574,388:392),1)
     ThetaT2_rec =
rcm2rgrid_Wrap(lat2(:,389:391),lon2(:,389:391),ThetaT2,LAT2,LON2,1)
     ThetaT3_rec =
rcm2rgrid_Wrap(lat3(:,389:391),lon3(:,389:391),ThetaT3,LAT3,LON3,1)

     levt_rec =
rcm2rgrid_Wrap(lat(:,389:391),lon(:,389:391),levt,lat(:,390),lon(574,389:391),1)

     dThdt = (ThetaT3_rec - ThetaT2_rec)/(2*36*60)
     dThDT = dThdt + US_rec*dThdx + VS_rec*dThdy + WS_rec(1:,:,:) * dThdz

     DBgrad = grad_latlon_cfd(dThDT, LAT, LON, False,False)

     dDBdy = DBgrad[0]
     dDBdx = DBgrad[1]
     DB = deltaTh *(dThdx * dDBdx + dThdy * dDBdy)
     DB_avg = dim_avg_n(DB,2)

All values calculated for DB_avg are missing values, even though deltaTh
*(dThdx * dDBdx) and deltaTh *(dThdy * dDBdy) alone do render reasonable
values for a lot of the points, but somehow when summed together,
everything goes wrong. In the case of Sh, Conv and Tilt, there are also a
significant number of missing values (more than I expected) but NCL can
plot the output data eventually. But that is not true for DB_avg.
I suspect the core problem should lie in converting curvilinear to
rectilinear grid, but I'm confused why summing two arrays should result in
missing values for every single dimension.

I'm also not sure if I'm calculating the gradients correctly. I really
appreciate your help.

Bests,
Melissa

-- 
*Melissa Kazemi Rad*
*Atmospheric Sciences Department*
*Rutgers University*
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20181218/18a7d8ef/attachment.html>