[ncl-talk] How to define date for WRF-chem to get average daily of a variable ?
Setareh Rahimi
setareh.rahimi at gmail.com
Sat Oct 8 07:05:51 MDT 2022
Dear Dennis,
I executed the script you send. It also worked for me. But may I ask you
please to execute the attached script?
Thanks in advance,
Best wishes,
On Fri, Oct 7, 2022 at 10:57 PM Dennis Shea <shea at ucar.edu> wrote:
> Yes, I did "execute the entire script" I sent you.
>
> %> /project/cas/shea/WRF> ncl setareh_chem.ncl >! setareh_chem.OUTPUT
> ====================================
>
> My script has no "x"
>
> I have no idea where or what the variable "x" is.
>
> Did you execute the script I sent you ?
>
>
>
> On Fri, Oct 7, 2022 at 1:02 PM Setareh Rahimi <setareh.rahimi at gmail.com>
> wrote:
>
>> Dear Dennis,
>> Did you execute the entire script?
>> Cause I still get errors! (attached)
>>
>> On Fri, Oct 7, 2022 at 9:31 PM Dennis Shea <shea at ucar.edu> wrote:
>>
>>> My script follows.. There were no problems. The output is attached.
>>>
>>> ++++++++++++++++++++++++++++++++++++++++++++++++++++
>>>
>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
>>>
>>> f = addfile("setareh_chem.nc","r")
>>> print(f) ; same as "nmcdump -h"
>>>
>>> Times = f->Times ; [Time | 41] x [DateStrLen |
>>> 19] (character)
>>>
>>> ; For fun: Examine "Times" in different formats
>>>
>>> Time_0 = wrf_times_c( Times, 0 ) ; "hours since" initial time on
>>> file (double)
>>> ; Time_0
>>> is recognized by cd_calendar
>>> ;Time_1 = wrf_times_c( Times, 1 ) ; "hours since 1901-01-01
>>> 00:00:00" (double)
>>> ;Time_2 = wrf_times_c( Times, 2 ) ; yyyymmddhhmnss
>>> (double)
>>> Time_3 = wrf_times_c( Times, 3 ) ; yyyymmddhh
>>> (integer)
>>>
>>> ;print(Times)
>>> ;print(Time_0)
>>> ;print(Time_1)
>>> ;print(Time_2)
>>> print(Time_3) ; easy on the eyes
>>> print("------------------------------------")
>>>
>>> a2 = f->TAUAER2 ; (Time, south_north, west_east)
>>> a2&Time = Time_0 ; associate the 'Time'
>>> coordinate with the variable using standard NCL & syntax
>>> printVarSummary(a2)
>>> print("------------------------------------")
>>>
>>> opt_cdv = True ; option for calculate_daily_values (cdv)
>>> opt_cdv at nval_crit = 4 ; default is 1
>>>
>>> a2Day = calculate_daily_values (a2, "avg", 0, opt_cdv)
>>> printVarSummary(a2Day)
>>> printMinMax (a2Day,1)
>>>
>>> On Fri, Oct 7, 2022 at 9:28 AM Setareh Rahimi <setareh.rahimi at gmail.com>
>>> wrote:
>>>
>>>> Dear Dennis,
>>>> I sent you the file via transfer.pcloud.com. I could not send it here
>>>> due to the large size of the file.
>>>> Please check the script with the file and let me know your idea.
>>>> Many thanks for your help.
>>>> Best wishes,
>>>>
>>>> On Fri, Oct 7, 2022 at 3:18 PM Dennis Shea <shea at ucar.edu> wrote:
>>>>
>>>>> I do not understand 'x' either
>>>>>
>>>>> "I wonder about variable x; I did not define such a variable."
>>>>>
>>>>> ncl 0> f = addfile(".....","r") ; I do not have your file
>>>>> ncl 1> print(f) ; like
>>>>> ncdump -h
>>>>> ncl 2> Times = f->Times ; Times(Time,
>>>>> DateStrLen) (type character)
>>>>> ncl 4> Time = wrf_times_c( Times, 0 ) ; "hours since" initial
>>>>> time on file (double); units recognized by cd_calendar
>>>>> ncl 5> print(Time)
>>>>>
>>>>> ncl 7> T = f->T
>>>>> ncl 8> printVarSummary(T)
>>>>>
>>>>> ncl 9> T&Time = Time ; associate time
>>>>> coordinate
>>>>> ncl 10> printVarSummary(T)
>>>>>
>>>>> Variable: T
>>>>> Type: float
>>>>> Total Size: 19683000 bytes
>>>>> 4920750 values
>>>>> Number of Dimensions: 4
>>>>> Dimensions and sizes: [Time | ??] x [bottom_top | 27] x [south_north |
>>>>> 81] x [west_east | 90]
>>>>> Coordinates:
>>>>> Time: [ 0.. ?]
>>>>> Number Of Attributes: 5
>>>>> FieldType : 104
>>>>> MemoryOrder : XYZ
>>>>> description : perturbation potential temperature (theta-t0)
>>>>> units : K
>>>>> stagger :
>>>>> ncl 11> Tavg = calculate_daily_values (T, "avg", 0, True)
>>>>> ncl 12> printVarSummary(Tavg)
>>>>>
>>>>> Variable: Tavg
>>>>> Type: float
>>>>> Total Size: 1574640 bytes
>>>>> 393660 values
>>>>> Number of Dimensions: 4
>>>>> Dimensions and sizes: [Time | ...] x [bottom_top | 27] x [south_north
>>>>> | 81] x [west_east | 90]
>>>>> Coordinates:
>>>>> Time: [ 0.. 12]
>>>>> Number Of Attributes: 8
>>>>> _FillValue : 9.96921e+36
>>>>> FieldType : 104
>>>>> MemoryOrder : XYZ
>>>>> description : perturbation potential temperature (theta-t0)
>>>>> units : K
>>>>> stagger :
>>>>> Time : 0
>>>>> NCL_tag : calculate_daily_values
>>>>>
>>>>> On Fri, Oct 7, 2022 at 5:16 AM Setareh Rahimi <
>>>>> setareh.rahimi at gmail.com> wrote:
>>>>>
>>>>>> Dear Dennis,
>>>>>> So many thanks for your advice. I made corrections to my script as
>>>>>> you suggested, and attempted to calculate and plot the average AOD for
>>>>>> 2011-06-28 but I still get errors. I attached the script and the errors.
>>>>>>
>>>>>> As you can see from the error message, NCL says :fatal: No
>>>>>> coordinate variable exists for dimension (Time) in variable (x)
>>>>>> I wonder about variable x; I did not define such a variable.
>>>>>>
>>>>>> Please advise me on how to sort this issue out.
>>>>>> Many thanks in advance,
>>>>>> Best wishes,
>>>>>>
>>>>>> On Fri, Oct 7, 2022 at 6:02 AM Dennis Shea <shea at ucar.edu> wrote:
>>>>>>
>>>>>>> As noted in the documentation, *calculate_daily_values*
>>>>>>> <https://www.ncl.ucar.edu/Document/Functions/Contributed/calculate_daily_values.shtml>
>>>>>>> requires that the variable have a **time coordinate** that is recognized by
>>>>>>> cd_calendar.
>>>>>>>
>>>>>>> Times = f->Times ; Times(Time, DateStrLen) (type character)
>>>>>>> Time = *wrf_times_c*( Times, 0 ) ; "hours since" initial time on file (double); units recognized by cd_calendar
>>>>>>>
>>>>>>> a2 = f->TAUAER ; (Time, south_north, west_east)
>>>>>>> a2&Time = Time ; associate the 'Time' coordinate with the
>>>>>>> variable using standard NCL & syntax
>>>>>>> printVarSummary(a2)
>>>>>>> print("------------------------------------")
>>>>>>>
>>>>>>> a2Day = calculate_daily_values (a2, "avg", 0, opt)
>>>>>>> printVarSummary(a2Day)
>>>>>>> printMinMax (a2Day,1)
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> On Wed, Oct 5, 2022 at 9:30 AM Setareh Rahimi via ncl-talk <
>>>>>>> ncl-talk at mailman.ucar.edu> wrote:
>>>>>>>
>>>>>>>> Dear NCL users,
>>>>>>>>
>>>>>>>> I am trying to plot the daily average of AOD (aerosol optical
>>>>>>>> depth) from WRF-chem outputs. However, I can not define the time so that
>>>>>>>> NCL plots the daily average of AOD.(the script has been attached, and the
>>>>>>>> model output is too large and I just used ncdump to show the header of some
>>>>>>>> varibles ).
>>>>>>>> I tried the following commands:
>>>>>>>>
>>>>>>>>
>>>>>>>> times = wrf_user_getvar(f,"times",-1)
>>>>>>>> printVarSummary(times)
>>>>>>>> print(times)
>>>>>>>> ymdh = cd_calendar(a2×, -2)
>>>>>>>> print(ymdh)
>>>>>>>>
>>>>>>>> but faced an error.
>>>>>>>> Would you please kindly advise me on how can I calculate the daily
>>>>>>>> average of AOD from the hourly output of WRF-chem?
>>>>>>>> Many thanks in advance,
>>>>>>>> Best wishes,
>>>>>>>> --
>>>>>>>> S.Rahimi
>>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>>> ncl-talk mailing list
>>>>>>>> ncl-talk at mailman.ucar.edu
>>>>>>>> List instructions, subscriber options, unsubscribe:
>>>>>>>> https://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>>>>>
>>>>>>>
>>>>>>
>>>>>> --
>>>>>> S.Rahimi
>>>>>>
>>>>>>
>>>>
>>>> --
>>>> S.Rahimi
>>>>
>>>>
>>
>> --
>> S.Rahimi
>>
>>
--
S.Rahimi
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20221008/c0abc5f5/attachment-0001.htm>
-------------- next part --------------
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/WRF_contributed.ncl"
;***************************************************************************************************************************
f = addfile ("wrfout_d01_2011-07-01_00:00:12", "r")
print(f) ; same as "nmcdump -h"
Times = f->Times ; [Time | 41] x [DateStrLen | 19] (character)
Time_0 = wrf_times_c( Times, 0 ) ; "hours since" initial time on file (double)
; Time_0 is recognized by cd_calendar
print("********************************************************************************************************************")
wks = gsn_open_wks("pdf" ,"WRF_aod_Daily-Average") ; ps,pdf,x11,ncgm,eps
gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map
res = True ; plot mods desired
res at cnConstFEnableFill = True
res at gsnMaximize = True
res at gsnSpreadColors = True
res at cnFillOn = True
res at cnLinesOn = False
res at cnLineLabelsOn = False
WRF_map_c(f, res, 0)
res at tfDoNDCOverlay = True
res at pmTickMarkDisplayMode = "Always"
;***************************************************************************************************************************
ymd = cd_calendar(Time_0, -2)
YMD = 20110701
mt = ind(YMD.eq.ymd)
a2 = f->TAUAER2(mt,:,:,:) ; (Time, south_north, west_east)
a2&Time = Time_0 ; associate the 'Time' coordinate with the variable using standard NCL & syntax
printVarSummary(a2)
print(ymd)
print("********************************************************************************************************************")
opt_cdv = True ; option for calculate_daily_values (cdv)
opt_cdv at nval_crit = 4 ; default is 1
a2Day = calculate_daily_values (a2, "avg", 0, opt_cdv)
printVarSummary(a2Day)
printMinMax (a2Day,1)
b3 = f->TAUAER3(mt,:,:,:)
b3Day = calculate_daily_values (b3, "avg", 0, opt_cdv)
printVarSummary(bDay)
printMinMax (bDay,1)
;********************************************************************************************************************
;********************************************************************************************************************
angstrom_exponent = -(log(a2Day)-log(b3Day))/log(400/600)
printVarSummary(angstrom_exponent)
AOD550_3D = a2Day * ((400/550)^angstrom_exponent)
printVarSummary(AOD550_3D)
AOD550_2D = dim_sum_n_Wrap(AOD550_3D,0)
if (any(isnan_ieee(AOD550_2D))) then
value = 9.96921e+36
replace_ieeenan (AOD550_2D,value,0)
AOD550_2D at _FillValue = value
end if
printVarSummary(AOD550_2D)
printMinMax (AOD550_2D,1)
res at tiMainString = "WRF-CHEM (aod_550) " + times(mt)
res at gsnLeftString = a at description
plot = gsn_csm_contour_map(wks,AOD550_2D(:,:),res)
end
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Screen Shot 1401-07-16 at 4.28.26 PM.png
Type: image/png
Size: 145175 bytes
Desc: not available
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20221008/c0abc5f5/attachment-0001.png>
More information about the ncl-talk
mailing list