<div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div class="gmail_default" style="font-family:georgia,serif">Hi,</div><div class="gmail_default" style="font-family:georgia,serif"><br></div><div class="gmail_default" style="font-family:georgia,serif">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!</div><div class="gmail_default" style="font-family:georgia,serif"><br></div><div class="gmail_default" style="font-family:georgia,serif">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.</div><div class="gmail_default" style="font-family:georgia,serif"><br></div><div class="gmail_default" style="font-family:georgia,serif"> begin<br><br>  a = addfile("wrfout_d02_2017-07-17_00:00:00","r")<br><br>  type = "pdf"<br><br>  wks = gsn_open_wks(type,"Kinematics")<br>  times  = wrf_user_list_times(a)<br>  time   = 40<br><br>  nplots     = 4<br>  plot       = new(nplots,"graphic")<br>   res                      = True<br>   res@tiMainOn             = False</div><div class="gmail_default" style="font-family:georgia,serif">;--------------------------- Shear -----------------------------------------;<br>     P   = wrf_user_getvar(a,"P",time)<br>     PB  = wrf_user_getvar(a,"PB",time)<br>     T   = wrf_user_getvar(a,"T",time)<br>     U   = wrf_user_getvar(a,"U",time)<br>     V   = wrf_user_getvar(a,"V",time)<br>     lon = wrf_user_getvar(a,"XLONG",time)<br>     lat = wrf_user_getvar(a,"XLAT",time)<br><br>     T = T + 300.<br>     P = P + PB<br><br>     TH = pot_temp(P, T, -1, False)<br>     Theta = TH(:,:,389:391)<br>     US = U(:,:,389:391)<br>     VS = V(:,1:,389:391)  ; The domain is 900*900 grid points.<br>     LAT = lat(:,390)<br>     LON = lon(574,389:391)<br><br>;------------- Convert Curvilinear to Rectilinear -----------<br>     Theta_rec = rcm2rgrid_Wrap(lat(:,389:391),lon(:,389:391),Theta,LAT,LON,1)<br>     US_rec    = rcm2rgrid_Wrap(lat(:,389:391),lon(:,389:391),US,LAT,LON,1)<br>     VS_rec    = rcm2rgrid_Wrap(lat(:,389:391),lon(:,389:391),VS,LAT,LON,1)<br>;-------------------- Gradients -----------------------------<br>     Thgrad = grad_latlon_cfd(Theta_rec,LAT,LON,False,False)<br>     Ugrad  = grad_latlon_cfd(US_rec,LAT,LON,False,False)<br>     Vgrad  = grad_latlon_cfd(VS_rec,LAT,LON,False,False)<br>     dThdx = Thgrad[1]<br>     dThdy = Thgrad[0]<br><br>     deltaTh = ((dThdx ^2) + (dThdy ^2))^(-0.5)<br>     dUdy  = Ugrad[0]<br>     dVdx  = Vgrad[1]<br><br>     Sh = - deltaTh *(dThdx * dThdy*(dUdy + dVdx))<br></div><div class="gmail_default" style="font-family:georgia,serif">    ; The above shear term calculated for a S-N cross section centered on i = 390<br></div><div class="gmail_default" style="font-family:georgia,serif">     Sh_avg = dim_avg_n(Sh,2)<br>;------------------Confluence ---------------------------------<br>     dVdy  = Vgrad[0]<br>     dUdx  = Ugrad[1]<br>     Con = - deltaTh *(((dThdy ^2) * dVdy) + ((dThdx ^2) * dUdx))<br>     Con_avg = dim_avg_n(Con,2)<br>;------------------- Tilting ------------------------<br>     W   = wrf_user_getvar(a,"W",time)<br>     Z   = wrf_user_getvar(a,"z",time)<br>     lev = Z(:,:,389:391)<br>     WS = W(:,:,389:391)<br>     WS_rec = rcm2rgrid_Wrap(lat(:,389:391),lon(:,389:391),WS,LAT,LON,1)<br>     lev_rec = rcm2rgrid_Wrap(lat(:,389:391),lon(:,389:391),lev,LAT,LON,1)<br><br>     dThdz = center_finite_diff_n(Theta_rec,lev_rec,False,0,0)<br><br>     Wgrad  = grad_latlon_cfd(WS_rec,LAT,LON,False,False)<br>     dWdy  = Wgrad[0]<br>     dWdx  = Wgrad[1]<br><br>     Tilt = - deltaTh *(dThdz * (dWdy(1:,:,:) * dThdy + dWdx(1:,:,:) * dThdx))<br>     Tilt_avg = dim_avg_n(Tilt,2)<br>;------------------------ Diabatic Heating Term ----------------------<br><br><br>     P2  = wrf_user_getvar(a,"P",time-1)<br>     P3  = wrf_user_getvar(a,"P",time+1)<br><br>     PB2  = wrf_user_getvar(a,"PB",time-1)<br>     PB3  = wrf_user_getvar(a,"PB",time+1)<br><br>     lat2 = wrf_user_getvar(a,"XLAT",time-1)<br>     lat3 = wrf_user_getvar(a,"XLAT",time+1)<br><br>     lon2 = wrf_user_getvar(a,"XLONG",time-1)<br>     lon3 = wrf_user_getvar(a,"XLONG",time+1)<br><br>     T2  = wrf_user_getvar(a,"T",time-1)<br>     T3  = wrf_user_getvar(a,"T",time+1)<br><br>     Tr2 = T2(:,:,389:391) + 300.<br>     Tr3 = T3(:,:,389:391) + 300.<br><br>     Pr2 = P2(:,:,389:391) + PB2(:,:,389:391)<br>     Pr3 = P3(:,:,389:391) + PB3(:,:,389:391)<br><br>     levt = Z(:,:,389:391)<br><br>     ThetaT2 = pot_temp(Pr2, Tr2, -1, False)<br>     ThetaT3 = pot_temp(Pr3, Tr3, -1, False)<br>     LAT2 = lat2(:,390)<br>     LON2 = lon2(574,389:391)<br><br>     LAT3 = lat3(:,390)<br>     LON3 = lon3(574,389:391)<br><br>rcm2rgrid_Wrap(lat(:,388:392),lon(:,388:392),ThetaT1,lat(:,390),lon(574,388:392),1)<br>     ThetaT2_rec = rcm2rgrid_Wrap(lat2(:,389:391),lon2(:,389:391),ThetaT2,LAT2,LON2,1)<br>     ThetaT3_rec = rcm2rgrid_Wrap(lat3(:,389:391),lon3(:,389:391),ThetaT3,LAT3,LON3,1)<br><br>     levt_rec = rcm2rgrid_Wrap(lat(:,389:391),lon(:,389:391),levt,lat(:,390),lon(574,389:391),1)<br><br>     dThdt = (ThetaT3_rec - ThetaT2_rec)/(2*36*60)</div><div class="gmail_default" style="font-family:georgia,serif">     dThDT = dThdt + US_rec*dThdx + VS_rec*dThdy + WS_rec(1:,:,:) * dThdz<br><br>     DBgrad = grad_latlon_cfd(dThDT, LAT, LON, False,False)<br><br>     dDBdy = DBgrad[0]<br>     dDBdx = DBgrad[1]<br>     DB = deltaTh *(dThdx * dDBdx + dThdy * dDBdy)</div><div class="gmail_default" style="font-family:georgia,serif">     DB_avg = dim_avg_n(DB,2)<br></div><div class="gmail_default" style="font-family:georgia,serif"><br></div><div class="gmail_default" style="font-family:georgia,serif">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.</div><div class="gmail_default" style="font-family:georgia,serif">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.</div><div class="gmail_default" style="font-family:georgia,serif"><br></div><div class="gmail_default" style="font-family:georgia,serif">I'm also not sure if I'm calculating the gradients correctly. I really appreciate your help.</div><div class="gmail_default" style="font-family:georgia,serif"><br></div><div class="gmail_default" style="font-family:georgia,serif">Bests,</div><div class="gmail_default" style="font-family:georgia,serif">Melissa<br></div><div class="gmail_default" style="font-family:georgia,serif"><br clear="all"></div><br>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><font face="georgia, serif"><i>Melissa Kazemi Rad</i></font><div><font face="georgia, serif"><i>Atmospheric Sciences Department</i></font></div><div><font face="georgia, serif"><i>Rutgers University</i></font></div></div></div></div></div></div></div></div></div></div></div></div></div></div></div></div></div></div></div>