[ncl-talk] Temp advection, wrf_contour fails
Dennis Shea
shea at ucar.edu
Thu Apr 21 10:24:51 MDT 2016
I think you would have to ask wrf-help at ucar.edu
I don't think this is an NCL-issue per-se. NCL is providing an interface to
fortran subroutines provided by the WRF project.
You state that the "latitude off by ~1 deg" ? Not sure how you get this.
I'm not a WRF user/expert so I can not be of much further help.
==
Yes ... people have calculated gradients on curvilinear grids but they use
ESMF regridding to interpolate to a rectilinear grid; perform the gradient
calculations and then reinterpolate back to the original curvilinear grid.
Not gradients related per se but it give you an idea of how to proceed ...
https://www.ncl.ucar.edu/Applications/ESMF.shtml
Example 5
Example 29
Example 30 is for a NARR grid (curvilinear). It illustrates that the errors
of back and forth interpolation are small.: original values are 200-300
.... differeines are 0.05. That's pretty good.
D
On Thu, Apr 21, 2016 at 9:55 AM, Ryan Connelly <rconne01 at gmail.com> wrote:
> Yeah that makes sense why wrf_user_intrp3d would fail. Your second
> suggestion got me lat/lon minima and maxima that are close, but not close
> enough, to the corners of my domain (latitude off by ~1 deg).
>
> I'm starting to think this isn't physically possible in NCL. The fact
> that WRF is a curvilinear grid means that I couldn't even manually code up
> conditionals and loops to calculate the gradient components, right? Since
> advancing from i=2 to i=4, e.g., isn't moving due west-east.
>
> Has no one ever computed advections of WRF variables in NCL before?
>
>
>
> On Wed, Apr 20, 2016 at 6:30 PM, Dennis Shea <shea at ucar.edu> wrote:
>
>> Likely, I was not clear. The following file is needed.
>>
>> wrfout_d01_2016-02-04_12:00:00
>>
>> ===
>> You should look at:
>> lon_plane = wrf_user_intrp3d(lon,p,"h",pressure,0.,False)
>> lat_plane = wrf_user_intrp3d(lat,p,"h",pressure,0.,False)
>>
>> printVarSummary(lon_plane)
>> printMinMax(lon_plane,0)
>>
>> printVarSummary(lat_plane)
>> printMinMax(lat_plane,0)
>>
>> Really, since lat/lon have no vertical extent .... how can you use
>> wrf_user_intrp3d?
>> https://www.ncl.ucar.edu/Document/Functions/WRF_arw/wrf_user_intrp3d.shtml
>>
>> The 1st arguments is:
>>
>> *var3d*
>>
>> Data on model levels that will be interpolated. The rightmost dimensions
>> of this array is *nz* x *ny* x *nx*.
>>
>> I am not a WRF user ... perhaps someone more knowledgeable can provide an
>> answer?
>>
>> Likely, you could 'trick' lat/lon to multiple dimensions via
>>
>> lat3d = conform(p, lat, (/1,2/))
>> lon3d = conform(p, lon, (/1,2/))
>>
>> lon_plane = wrf_user_intrp3d(lon3d,p,"h",pressure,0.,False)
>> lat_plane = wrf_user_intrp3d(lat3d,p,"h",pressure,0.,False)
>>
>>
>>
>>
>> ===
>> re: "realistic values" for dTdx instead of inf's:
>>
>> (0) min=-15.1425 max=18.5633
>> (0) min=-24.1597 max=15.1911
>>
>> ===
>> The actual gradients will likely *not* be these nice numbers.
>>
>> The (say): (T(n+1)-T(n-1))/(LAT(n+1)-LAT(n-1))
>>
>> (20-10)/(15-12.5)
>>
>> would yield 10/2.5 = 4 [degC/degLat] or degF in your case
>>
>> however, most gradients are (say) degC/m_or_km
>>
>> so that 2.5 degrees of lat would have to be converted to m or km.
>>
>> ===
>> You are aware that the WRF grid is curvilinear ... not rectilinear?
>>
>>
>>
>> On Wed, Apr 20, 2016 at 11:49 AM, Ryan Connelly <rconne01 at gmail.com>
>> wrote:
>>
>>> Hi Dennis,
>>>
>>> My file is now in ../incoming/wrf_temp_adv_rconne01.ncl.
>>>
>>> It doesn't like dimension reduction because the resulting argument isn't
>>> the right size:
>>>
>>> fatal:center_finite_diff_n: r must either be a scalar, a 1D array the
>>> same length as the dim-th dimemsion of q, or the same size as q
>>>
>>>
>>> It also throws a fatal on the declaration of lon_plane and/or lat_plane:
>>>
>>> fatal:Subscript out of range, error in subscript #0
>>> fatal:An error occurred reading dims
>>> fatal:["Execute.c":8128]:Execute: Error occurred at or near line 224 in
>>> file $NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl
>>>
>>> fatal:["Execute.c":8128]:Execute: Error occurred at or near line 97 in
>>> file wrf_temp_adv.ncl
>>>
>>>
>>> My best success so far has been with your first suggestion, merely
>>> changing the last arg of the dTDx declaration from 0 to 1, which gave me
>>> realistic values for dTdx instead of inf's:
>>>
>>> (0) min=-15.1425 max=18.5633
>>> (0) min=-24.1597 max=15.1911
>>>
>>> But that ultimately still hangs up, though interestingly enough, it now
>>> gets one step farther. Before, wrf_contour was where it stopped. Now, it
>>> gets past wrf_contour, but fails on plot=.
>>>
>>> Thanks,
>>> Ryan
>>>
>>> On Wed, Apr 20, 2016 at 1:19 PM, Ryan Connelly <rconne01 at gmail.com>
>>> wrote:
>>>
>>>> Hi Mary,
>>>>
>>>> That shows that dTdx is the problem. (dTdx is first, dTdy second:)
>>>>
>>>> (0) min=-inf max=inf
>>>> (0) min=-24.1597 max=15.1911
>>>>
>>>> I'll work on Dennis's suggestion now.
>>>>
>>>> Ryan
>>>>
>>>> On Wed, Apr 20, 2016 at 11:59 AM, Dennis Shea <shea at ucar.edu> wrote:
>>>>
>>>>> NOTE: I have not looked at this in detail.
>>>>>
>>>>> [1]
>>>>> You could not use 'advect_variable' because it requires "[t]he array *must
>>>>> be global and ordered south to north*." Also, this assumes a
>>>>> rectilinear grid (eg: fixed or gaussian)
>>>>>
>>>>> Your data are not global :-(
>>>>> Yourhorizontal grid is not rectilinear.
>>>>>
>>>>> Note: I could change ''advect_variable" which uses spherical harmonics
>>>>> (hence, the global grid requirement) to derive the horizontal and
>>>>> latitudinal gradients tusing the 6.4.0
>>>>>
>>>>> https://www.ncl.ucar.edu/Document/Functions/Contributed/grad_latlon_cfd.shtml
>>>>> BUT that would take a day or two.
>>>>>
>>>>> [2]
>>>>> center_finite_diff "the issue" .. yes and no
>>>>>
>>>>> The error occurs in center_finite_diff but the issue is, I speculate,
>>>>> the arguments. Specifically, the use of lat and lon
>>>>>
>>>>> ====
>>>>> lon = wrf_user_getvar(a,"lon",it)
>>>>> lat = wrf_user_getvar(a,"lat",it)
>>>>>
>>>>> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>>>>> ; Interpolate to isobaric planes
>>>>>
>>>>> pressure = 850.
>>>>>
>>>>> tc_plane = wrf_user_intrp3d(tc,p,"h",pressure,0.,False)
>>>>> z_plane = wrf_user_intrp3d(z_dec,p,"h",pressure,0.,False)
>>>>> u_plane = wrf_user_intrp3d( u,p,"h",pressure,0.,False)
>>>>> v_plane = wrf_user_intrp3d( v,p,"h",pressure,0.,False)
>>>>>
>>>>>
>>>>> ; Define variables needed to calculate temp advection
>>>>>
>>>>> dTdx = center_finite_diff_n(tc_plane,lon,False,0,0)
>>>>> dTdy = center_finite_diff_n(tc_plane,lat,False,0,0)
>>>>> ===
>>>>>
>>>>> Maybe,
>>>>>
>>>>>
>>>>> dTdx = center_finite_diff_n(tc_plane,lon,False,0,1) ; <==
>>>>> dTdy = center_finite_diff_n(tc_plane,lat,False,0,0)
>>>>>
>>>>> Or [281] x [353]=>nlat,mlon
>>>>>
>>>>> dTdx =
>>>>> center_finite_diff_n(tc_plane,lon(nlat/2,:),False,0,0)
>>>>> dTdy = center_finite_diff_n(tc_plane,lat(:,mlon/2),False,0,0)
>>>>>
>>>>> Note that NCL's 'dimension reduction' (elimination of degenerate
>>>>> dimensions), lon(nlat/2,:) and lat(:,mlon/2) are one dimensional.
>>>>> Use printVarSummary( lon(nlat/2,:) ) to verify
>>>>>
>>>>> [3]
>>>>>
>>>>> [2] is a thumb in the dike approach. I think it would be best (if
>>>>> possible) to calculate
>>>>>
>>>>> lon_plane =
>>>>> lat_plane =
>>>>>
>>>>> dTdx = center_finite_diff_n(tc_plane,lon_plane,False,0,0)
>>>>> dTdy = center_finite_diff_n(tc_plane,lat_plane,False,0,0)
>>>>>
>>>>> [4]
>>>>> Maybe you could send a sample filw (one only) to
>>>>>
>>>>> ftp ftp.cgd.ucar.edu
>>>>> anonymous
>>>>> your_email
>>>>> cd incoming
>>>>> put .....
>>>>> quit
>>>>>
>>>>> then let ncl_talk know the name of the file.
>>>>>
>>>>> On Tue, Apr 19, 2016 at 4:21 PM, Mary Haley <haley at ucar.edu> wrote:
>>>>>
>>>>>> Hi Ryan,
>>>>>>
>>>>>> That makes sense why you can't upgrade.
>>>>>>
>>>>>> As for the "inf", this is definitely an issue. To further determine
>>>>>> where the problem is, I suggest doing a printMinMax on the individual
>>>>>> variables in the calculation:
>>>>>>
>>>>>> printMinMax(dTdx,0)
>>>>>> printMinMax(dTdy,0)
>>>>>>
>>>>>> etc.
>>>>>>
>>>>>> You may want to do this on every variable that you are getting with
>>>>>> "wrf_user_getvar" to make sure those all look okay.
>>>>>>
>>>>>> If it looks like center_finite_diff is the issue, then revisit the
>>>>>> documentation to make sure all the input parameters are being input
>>>>>> correctly:
>>>>>>
>>>>>>
>>>>>> http://www.ncl.ucar.edu/Document/Functions/Built-in/center_finite_diff.shtml
>>>>>>
>>>>>> --Mary
>>>>>>
>>>>>> On Tue, Apr 19, 2016 at 3:53 PM, Ryan Connelly <rconne01 at gmail.com>
>>>>>> wrote:
>>>>>>
>>>>>>> Hi Mary,
>>>>>>>
>>>>>>> I use NCL that's installed on a shared server where WRF is
>>>>>>> configured. I don't suspect I have permissions to remove an old version.
>>>>>>> I recently updated to 6.3.0 on my local machine, but then I have to copy
>>>>>>> the wrfout over, which is a bit time-consuming since I'm running a triple
>>>>>>> nest down to 1.333 km! :O
>>>>>>>
>>>>>>> These extra commands that you suggested certainly give a clue that
>>>>>>> I'm not calculating things right:
>>>>>>>
>>>>>>> Variable: temp_adv
>>>>>>> Type: float
>>>>>>> Total Size: 396772 bytes
>>>>>>> 99193 values
>>>>>>> Number of Dimensions: 2
>>>>>>> Dimensions and sizes: [281] x [353]
>>>>>>> Coordinates:
>>>>>>> Number Of Attributes: 1
>>>>>>> _FillValue : 9.96921e+36
>>>>>>> (0) min=-inf max=inf
>>>>>>>
>>>>>>>
>>>>>>> 281x353 is my grid size, so that's good, but the infinities probably
>>>>>>> are not...
>>>>>>>
>>>>>>> Thanks,
>>>>>>> Ryan
>>>>>>>
>>>>>>> On Tue, Apr 19, 2016 at 10:38 AM, Mary Haley <haley at ucar.edu> wrote:
>>>>>>>
>>>>>>>> Hi Ryan,
>>>>>>>>
>>>>>>>> Is there a reason why you don't want to upgrade from V6.1.2? The
>>>>>>>> version you have is over 3 years old.
>>>>>>>>
>>>>>>>> I can't be sure of why your temp_adv variable is not plotting,
>>>>>>>> without actually being able to run the script.
>>>>>>>>
>>>>>>>> However, rather than printing the whole temp_adv array, what does
>>>>>>>> the following report:
>>>>>>>>
>>>>>>>> printVarSummary(temp_adv)
>>>>>>>> printMinMax(temp_adv,0)
>>>>>>>>
>>>>>>>> Sometimes this will give you a clue where the problem is, for
>>>>>>>> example, your min/max values are off scale, or the size of the array is not
>>>>>>>> what you were expecting.
>>>>>>>>
>>>>>>>> You can also try plotting temp_adv with a basic contour plot call:
>>>>>>>>
>>>>>>>> plot = gsn_csm_contour(wks,temp_adv,False)
>>>>>>>>
>>>>>>>> just to make sure that the problem isn't with wrf_contour itself.
>>>>>>>>
>>>>>>>> --Mary
>>>>>>>>
>>>>>>>>
>>>>>>>> On Mon, Apr 18, 2016 at 12:49 PM, Ryan Connelly <rconne01 at gmail.com
>>>>>>>> > wrote:
>>>>>>>>
>>>>>>>>> Hi,
>>>>>>>>>
>>>>>>>>> Running 6.1.2, so can't use advect_variable without upgrading. So
>>>>>>>>> instead I have...
>>>>>>>>>
>>>>>>>>> ;do it = 12,30,3 ; TIME LOOP
>>>>>>>>> do it = 12,13,1
>>>>>>>>>
>>>>>>>>> print("Working on time: " + times(it) )
>>>>>>>>> res at TimeLabel = times(it) ; Set Valid time to use on plots
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>>>>>>>>> ; First get the variables we will need
>>>>>>>>>
>>>>>>>>> tc = wrf_user_getvar(a,"tc",it) ; 3D tc
>>>>>>>>> td = wrf_user_getvar(a,"td",it) ; 3D td
>>>>>>>>> u = wrf_user_getvar(a,"ua",it) ; 3D U at mass points
>>>>>>>>> v = wrf_user_getvar(a,"va",it) ; 3D V at mass points
>>>>>>>>> p = wrf_user_getvar(a, "pressure",it) ; pressure is our
>>>>>>>>> vertical coordinate
>>>>>>>>> z = wrf_user_getvar(a,"z",it) ; Full model height in
>>>>>>>>> meters
>>>>>>>>> z_dec = z/10. ; Height in decameters
>>>>>>>>> lon = wrf_user_getvar(a,"lon",it)
>>>>>>>>> lat = wrf_user_getvar(a,"lat",it)
>>>>>>>>>
>>>>>>>>> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>>>>>>>>> ; Interpolate to isobaric planes
>>>>>>>>>
>>>>>>>>> pressure = 850.
>>>>>>>>>
>>>>>>>>> tc_plane = wrf_user_intrp3d(tc,p,"h",pressure,0.,False)
>>>>>>>>> z_plane = wrf_user_intrp3d(z_dec,p,"h",pressure,0.,False)
>>>>>>>>> u_plane = wrf_user_intrp3d( u,p,"h",pressure,0.,False)
>>>>>>>>> v_plane = wrf_user_intrp3d( v,p,"h",pressure,0.,False)
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> ; Define variables needed to calculate temp advection
>>>>>>>>>
>>>>>>>>> dTdx = center_finite_diff_n(tc_plane,lon,False,0,0)
>>>>>>>>> dTdy = center_finite_diff_n(tc_plane,lat,False,0,0)
>>>>>>>>>
>>>>>>>>> temp_adv = u_plane*dTdx + v_plane*dTdy
>>>>>>>>>
>>>>>>>>> print(temp_adv)
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>>>>>>>>>
>>>>>>>>> ; Plotting options for Tc
>>>>>>>>> opts = res
>>>>>>>>> print("Defined ops")
>>>>>>>>> opts at cnFillOn = True
>>>>>>>>> print("Turned fill on")
>>>>>>>>> opts at cnLabelMasking = True
>>>>>>>>> opts at cnInfoLabelOn = False
>>>>>>>>> opts at cnLineLabelPerimOn = False
>>>>>>>>> opts at ContourParameters = (/ -60., 60., 3./)
>>>>>>>>> print("Accepted contour parameters")
>>>>>>>>> ;opts at gsnSpreadColorEnd = -3 ; End third from the last
>>>>>>>>> color in color map
>>>>>>>>> contour_tc = wrf_contour(a,wks,temp_adv,opts) ; <- Breaks
>>>>>>>>> right here
>>>>>>>>> print("wrf_contour called")
>>>>>>>>> delete(opts)
>>>>>>>>> print("opts deleted")
>>>>>>>>>
>>>>>>>>> print("Got past plotting options for Tc")
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> ; MAKE PLOTS
>>>>>>>>> ;plot =
>>>>>>>>> wrf_map_overlays(a,wks,(/contour_tc,contour_td,contour_z/),pltres,mpres)
>>>>>>>>> plot = wrf_map_overlays(a,wks,(/contour_tc/),pltres,mpres)
>>>>>>>>> ;plot =
>>>>>>>>> wrf_map_overlays(a,wks,(/contour_td,vector/),pltres,mpres)
>>>>>>>>>
>>>>>>>>> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>>>>>>>>>
>>>>>>>>> end do ; END OF TIME LOOP
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> When it prints out the full grid of temp_adv, or random indexes of
>>>>>>>>> the array, I get reasonable values and no reason to suspect the grid is
>>>>>>>>> undefined anywhere. Yet as you can see from my print statements, the call
>>>>>>>>> to wrf_contour fails. It does not throw an error; it just hangs up. I
>>>>>>>>> checked with just tc_plane instead, and that completed fine, so it has to
>>>>>>>>> be an issue with the temp_adv grid I calculated.
>>>>>>>>>
>>>>>>>>> Thanks,
>>>>>>>>> Ryan
>>>>>>>>>
>>>>>>>>> --
>>>>>>>>> Ryan Connelly
>>>>>>>>> M.S. Student in Atmospheric Sciences, Stony Brook University
>>>>>>>>> B.S. in Meteorology with Minors in Mathematics and GIS, Valparaiso
>>>>>>>>> University
>>>>>>>>> rconne01 at gmail.com
>>>>>>>>> ryan.connelly at stonybrook.edu
>>>>>>>>>
>>>>>>>>> _______________________________________________
>>>>>>>>> ncl-talk mailing list
>>>>>>>>> ncl-talk at ucar.edu
>>>>>>>>> List instructions, subscriber options, unsubscribe:
>>>>>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>>>>>>
>>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> --
>>>>>>> Ryan Connelly
>>>>>>> M.S. Student in Atmospheric Sciences, Stony Brook University
>>>>>>> B.S. in Meteorology with Minors in Mathematics and GIS, Valparaiso
>>>>>>> University
>>>>>>> rconne01 at gmail.com
>>>>>>> ryan.connelly at stonybrook.edu
>>>>>>>
>>>>>>
>>>>>>
>>>>>> _______________________________________________
>>>>>> ncl-talk mailing list
>>>>>> ncl-talk at ucar.edu
>>>>>> List instructions, subscriber options, unsubscribe:
>>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>>>
>>>>>>
>>>>>
>>>>
>>>>
>>>> --
>>>> Ryan Connelly
>>>> M.S. Student in Atmospheric Sciences, Stony Brook University
>>>> B.S. in Meteorology with Minors in Mathematics and GIS, Valparaiso
>>>> University
>>>> rconne01 at gmail.com
>>>> ryan.connelly at stonybrook.edu
>>>>
>>>
>>>
>>>
>>> --
>>> Ryan Connelly
>>> M.S. Student in Atmospheric Sciences, Stony Brook University
>>> B.S. in Meteorology with Minors in Mathematics and GIS, Valparaiso
>>> University
>>> rconne01 at gmail.com
>>> ryan.connelly at stonybrook.edu
>>>
>>
>>
>
>
> --
> Ryan Connelly
> M.S. Student in Atmospheric Sciences, Stony Brook University
> B.S. in Meteorology with Minors in Mathematics and GIS, Valparaiso
> University
> rconne01 at gmail.com
> ryan.connelly at stonybrook.edu
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160421/e97897eb/attachment.html
More information about the ncl-talk
mailing list