[gradsusr] Total Snow Accumulation with GFS .25 Degree

Jeff Chabot jsc219 at gmail.com
Tue Jan 4 23:09:53 EST 2022 

Hello GrADS Users,

I solved the total snow accumulation equation for GFS .25 Deg.  Here is
what I came up and appears to work well:

TIME_STEPS = 81
'define precip = apcpsfc/25.4'
'define snow = precip * csnowsfc'

while ( t <= TIME_STEPS )
'set t 't
if (t=1)
'define rain = precip'
'define apcptot = precip'
'define snowtot = snow'
else
 'define rain = precip - precip(t-1)'
 'define apcptot = sum(precip,t=1,t-0) - sum(precip,t-1,t-1)'
 'define snowtot = (sum((precip(t-0) - precip(t-1))*csnowsfc,t=1,t-0,1)) *
10'
endif

'd rain'
'clear'
'd apcptot'
'clear'
'd snowtot'
'clear'

t = t+1
endwhile

Hopefully this helps others!  Jeff D, thanks for all your input!  If I can
find a single GRIB2 file that covers the forecast period, I will look into
switching methods, but for now the sdfopen method is working for me.

Sincerely,

Jeff C



On Mon, Jan 3, 2022 at 4:48 PM Jeff Chabot <jsc219 at gmail.com> wrote:

> Looking at that link, it appears they have one file per forecast hour as
> seen here:
>
>
> https://nomads.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/gfs.20220103/12/atmos/
>
> gfs.t12z.pgrb2.0p25.f000 <https://nomads.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/gfs.20220103/12/atmos/gfs.t12z.pgrb2.0p25.f000>
>
> Does that mean that I would need to download each forecast file that I want use?  That would require a lot of open commands in GrADS. Do they not store the whole GRIB data of GFS 0.25 Deg in one large file?
>
> Sincerely,
>
> Jeff C
>
>
>
>
>
>
>
> On Mon, Jan 3, 2022 at 12:22 PM Jeff Duda <jeffduda319 at gmail.com> wrote:
>
>> The NCEP nomads server for live GFS data, the directory of GFS data, is
>> at https://nomads.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/. They're
>> usually there for about 2 days before being purged.
>>
>> Yes, the files are somewhat large, so you may find delays with the
>> downloading process, but at least then you'd have the files locally.
>>
>> Yes, all your scripting would work identically with the GRIB formatted
>> data. As far as querying GRIB data outside of GrADS, either wgrib or wgrib2
>> are your best friends, and grib2ctl.pl and g2ctl can be used to write
>> control files from GRIB data. The only thing you need to pay attention to
>> is the GRIB version - there is a version 1 and a version 2. Version 2 is
>> more compressed than version 1, and it seems that most agencies have
>> replaced version 1 usage with version 2 over the past several years. wgrib2
>> and g2ctl are the particular tools that are used on GRIB2 data.
>>
>> I'm not sure about the grib filter data. It sounds like a promising
>> possibility, but I have not been able to get much out of it myself.
>>
>> Jeff Duda
>>
>> On Mon, Jan 3, 2022 at 9:45 AM Jeff Chabot <jsc219 at gmail.com> wrote:
>>
>>> I use and have been using sdfopen on the links from NOMADS to run all my
>>> scripts such as this link sinc 2003:
>>>
>>> https://nomads.ncep.noaa.gov/dods/gfs_0p25/gfs20220103/gfs_0p25_00z
>>>
>>> I have always done it that way and have been able to produce ptype,
>>> temp, and cloud products for multiple regions.  Just this year, I started
>>> downloading the files using sdfopen and fwrite so I would have the files
>>> locally and it makes for quicker runs and tests.  After speeding up the
>>> process, I thought it would be nice to try to create snow accumulation
>>> products.  I already have that working with NAM, GFS Ensembles, and HRRR.
>>> GFS .025 Deg however is just not cooperating.
>>>
>>> My latest attempt to get total snow accum to work with GFS .25 Deg looks
>>> like this:
>>>
>>> 'define precip = apcpsfc/25.4'
>>> 'define snow = precip * csnowsfc'
>>>
>>> while (1)
>>> 'set t 't
>>> if (t=1)
>>> 'define rain = precip'
>>> 'define apcptot = precip'
>>> 'define snowtot = snow'
>>> else
>>>  'define rain = precip - precip(t-1)'
>>>  'define apcptot = sum(precip,t=1,t-0) - sum(precip,t-1,t-1)'
>>>  fhr = (t-1) * 3
>>>  mod = math_fmod(fhr,6)
>>>  if (mod = 0)
>>> *  'define apcptot = sum(precip,t=1,t-0,2)' *doesn't work, makes precip
>>> too hot.
>>>   'define snowtot = (sum(snow(t-0) - snow(t-1),t=1,t-0,1)) * 10'
>>>  else
>>>   if (mod = 3)
>>> *   'define apcptot = sum(precip,t=1,t-1,2) + precip' *doesn't work,
>>> makes precip too hot.
>>>    'define snowtot = (sum(snow(t-0) - snow(t-1),t=1,t-0,1)) * 10'
>>>   endif
>>>  endif
>>> endif
>>>
>>> *bunch of code to make things look pretty*
>>> 'display rain'
>>>
>>> *bunch of code to make things look pretty*
>>> 'display apcptot'
>>>
>>> *bunch of code to make things look pretty*
>>> 'display snowtot'
>>>
>>> t = t+1
>>> if t > TIME_STEPS;
>>>   break;
>>> endif
>>> endwhile
>>>
>>> The above code results in negative snow at some "t" times, so I know
>>> it's wrong but I think it is close.  Using your sample equations results in
>>> overdone snow accum.  So, I have been tweaking the equations to find the
>>> correct solution.
>>>
>>> Examples of my testing can be found here (12Z 1/1/2022):
>>>
>>> https://jeffsweatherservice.com/grads/gfs_hd/test/snow_US_2.png
>>> https://jeffsweatherservice.com/grads/gfs_hd/test/precip_US_2.png
>>> https://jeffsweatherservice.com/grads/gfs_hd/test/ptype_US_2.png
>>>
>>> Raw grid analysts are here:
>>> https://jeffsweatherservice.com/grads/gfs_hd/test/csnowsfc_US_2.png
>>> https://jeffsweatherservice.com/grads/gfs_hd/test/apcpsfc_US_2.png
>>>
>>> Changing the "2" to any number up to "9" will show you the results I am
>>> seeing.
>>>
>>> As for GRIB files, I have never worked with them before.  Where would I
>>> obtain the files, would my existing scripts in GrADS work with GRIB files?
>>> Would I use the same open command in GrADS?  Can I download just the
>>> parameters that I want to use in the GRIB file?  I see the grib filter
>>> option on the NOMADS server, but I can't figure out what it does or how it
>>> works.
>>>
>>> Sincerely,
>>>
>>> Jeff C
>>>
>>>
>>> On Mon, Jan 3, 2022 at 12:38 AM Jeff Duda <jeffduda319 at gmail.com> wrote:
>>>
>>>> Jeff,
>>>> 'gxout fwrite' writes to a standard binary file, not a GRIB file. You
>>>> will not be able to use the wgrib products on such a file. Is there a
>>>> reason you need to write the snow data to a separate file? You can perform
>>>> math on arrays from a GRIB2 file. See this helpful write-up on wgrib2
>>>> tricks for inspiration:
>>>> https://ftp.cpc.ncep.noaa.gov/wd51we/wgrib2/tricks.wgrib2
>>>>
>>>> Jeff Duda
>>>>
>>>> On Thu, Dec 30, 2021 at 2:48 PM Jeff Chabot <jsc219 at gmail.com> wrote:
>>>>
>>>>> Hi Jeff,
>>>>>
>>>>> I have been manipulating your pseudo-code to try to get it to work and
>>>>> I am getting much closer. I will try this next: fhr = (t - 1)*3 in my
>>>>> equations.  I was actually trying (t+1)*3 and that didn't really help.  I
>>>>> also tried changing the start time in the sum line to t=2 (instead of t=1)
>>>>> to get past issues of apcpsfc at t=1.
>>>>>
>>>>> I am using the NOMADS server for the data, for example:
>>>>>
>>>>> https://nomads.ncep.noaa.gov/dods/gfs_0p25/gfs20211230/gfs_0p25_06z
>>>>>
>>>>> I use fwrite, display variables, and create a dat file.
>>>>>
>>>>> For example:
>>>>> ttime        = 1
>>>>> TIME_STEPS   = 129
>>>>> 'set fwrite gfs_hd'.dat'
>>>>> 'set gxout fwrite'
>>>>> while(1)
>>>>> 'set t 'ttime''
>>>>> 'set x 1 1440'
>>>>> 'set y 1 721'
>>>>> 'set lev 1000'
>>>>> 'd const(apcpsfc,0,-u)'
>>>>> 'd csnowsfc'
>>>>> ttime = ttime + 1
>>>>>
>>>>> if ttime > TIME_STEPS ;
>>>>>   break;
>>>>> endif
>>>>> endwhile
>>>>> 'disable fwrite'
>>>>>
>>>>> My ctl file has an "hr" in it:
>>>>> TDEF 129 linear 06Z30DEC2021 3hr
>>>>>
>>>>> The above method works great, producing other products.  I download
>>>>> other variables and plot them just fine.
>>>>>
>>>>> But, wgrib2 and g2ctl don't work with my grads generated dat file.  I
>>>>> tried both scripts, wgrib2 doesn't give me any output and g2ctl says it is
>>>>> not a girb file.
>>>>>
>>>>> Thanks again,
>>>>>
>>>>> Jeff C
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> On Thu, Dec 30, 2021 at 12:48 PM Jeff Duda <jeffduda319 at gmail.com>
>>>>> wrote:
>>>>>
>>>>>> Jeff,
>>>>>> My previous response included some pseudo-code...it was not meant to
>>>>>> be used verbatim. You will need to make adjustments for your specific
>>>>>> application.
>>>>>>
>>>>>> What is your source for the GFS data? Can you post output from wgrib2
>>>>>> -v or g2ctl ?
>>>>>>
>>>>>> As far as I've ever known, the time index in Grads always starts at
>>>>>> 1, regardless of the actual forecast length. The time axis is defined in
>>>>>> the control file...in the TDEF entry. So in reality, a formula converting
>>>>>> the time index in Grads to the valid forecast hour is fhr = (t - 1)*3
>>>>>> instead of what I previously said...unless of course your TDEF entry did
>>>>>> not include the 00-hr output from GFS.
>>>>>>
>>>>>> Remember, a modulus in math is simply a remainder. So (3 mod 6) = 3
>>>>>> even though the value of 3/6 is 0.5.
>>>>>>
>>>>>> This method does not depend on the initialization time.
>>>>>>
>>>>>> Jeff Duda
>>>>>>
>>>>>> On Thu, Dec 30, 2021 at 12:36 AM Jeff Chabot <jsc219 at gmail.com>
>>>>>> wrote:
>>>>>>
>>>>>>> Hi Jeff,
>>>>>>>
>>>>>>> The code below looks good and I thought it would solve my issue, but
>>>>>>> it didn't quite work.  Does it assume that t=1 is at the initial model run
>>>>>>> time such as 0Z, 6Z, 12Z, and 18Z?  Also, I had issues at t=1 with
>>>>>>> apcpsfc.  I had a lot of values there at: -39330708 based on a grid check.
>>>>>>> I believe that I corrected the issue by using const(apcpsfc,0,-u) when I
>>>>>>> download the .dat file to set all of those large negative numbers to 0.
>>>>>>>
>>>>>>> I am finding that I still need the following if statement or t=1
>>>>>>> will error out and t=2 will show no data:
>>>>>>>
>>>>>>> if (t=1|t=2)
>>>>>>>  'define snow = apcpsfc * csnowsfc'
>>>>>>> else
>>>>>>> ...
>>>>>>>
>>>>>>> Can you tell me if I am using the math_fmod properly?
>>>>>>>
>>>>>>> For example, let's say I'm using the 18Z model run, t=1 (18Z), t=2
>>>>>>> (21Z) and t=3 (0Z), etc...  At t = 3, then fhr = 3 * 3 = 9 and 9/6 is 1
>>>>>>> with a remainder of 3 so skip to the if (mod = 3).  Am I using that
>>>>>>> function correctly?  So, t =1 (say at 18Z), would be 3/6 = 0.5 which is not
>>>>>>> defined unless I use my if statement above.
>>>>>>>
>>>>>>> Thanks again,
>>>>>>>
>>>>>>> Jeff C
>>>>>>>
>>>>>>>
>>>>>>> On Wed, Dec 29, 2021 at 2:03 PM Jeff Duda <jeffduda319 at gmail.com>
>>>>>>> wrote:
>>>>>>>
>>>>>>>> csnowsfc represents a snow-mask only for falling precip in that
>>>>>>>> time step. It is not a fixed mask. Therefore, you will not be able to used
>>>>>>>> your proposed explicit formula for computing fallen snow at each time step.
>>>>>>>> You will need to sum the time-step snowfall to get run-total snowfall
>>>>>>>> accumulation using something like the following:
>>>>>>>>
>>>>>>>> fhr = t * 3
>>>>>>>> mod = math_fmod(fhr,6)
>>>>>>>> if (mod = 0)
>>>>>>>>  'define snow = sum(apcpsfc*csnowsfc,t=1,t-0,2)'
>>>>>>>> else
>>>>>>>>  if (mod = 3)
>>>>>>>>   'define snow = sum(apcpsfc*csnowsfc,t=1,t-1,2) + apcpsfc*csnowsfc'
>>>>>>>>  endif
>>>>>>>> endif
>>>>>>>>
>>>>>>>> 'd 10*snow/25.4'
>>>>>>>>
>>>>>>>> the ",2)" in the sum command refers to skipping one time slice
>>>>>>>> each, to avoid double counting those 3-hour increments. The second part of
>>>>>>>> the if-statement is to add the final 3-hour accumulated snowfall.
>>>>>>>>
>>>>>>>> There is a different way to do this, since those data files also
>>>>>>>> contain snowdepth as snodsfc. You could simply display 39.97*(snowdsfc -
>>>>>>>> snowdsfc(t-1)) as the difference in snow depth between 3-hour outputs
>>>>>>>> (scaled to inches from meters). Yes, there will be some negatives due to
>>>>>>>> snow melt, but you can easily mask those out. I would imagine this method
>>>>>>>> would be quite similar to the above method, as it is unlikely that there
>>>>>>>> would be much snow melt or removal at a grid point during which snow was
>>>>>>>> actively falling during a 3-hour window. And this latter formula is
>>>>>>>> certainly simpler and less computationally intensive.
>>>>>>>>
>>>>>>>> Jeff Duda
>>>>>>>>
>>>>>>>> On Tue, Dec 28, 2021 at 10:03 PM Jeff Chabot <jsc219 at gmail.com>
>>>>>>>> wrote:
>>>>>>>>
>>>>>>>>> Hello GrADS Users,
>>>>>>>>>
>>>>>>>>> I have been struggling to plot total snow accumulations with GFS
>>>>>>>>> .25 degree for the past week or so.  I know GFS does the following with
>>>>>>>>> precip accum: 0-3hr;0-6hr,6-9hr; 6-12hr, etc and I am trying to workaround
>>>>>>>>> that alternating precip totals.
>>>>>>>>>
>>>>>>>>> I was able to get this equation to work for total liquid precip,
>>>>>>>>> but I can't get it to work with snow:
>>>>>>>>>
>>>>>>>>> if t=1|t=2
>>>>>>>>> 'display 'precip'
>>>>>>>>> else
>>>>>>>>> 'display sum(precip,t-1,t+0) - sum(precip,t-1,t-1)'
>>>>>>>>> endif
>>>>>>>>>
>>>>>>>>> Where:
>>>>>>>>> precip = apcpsfc/25.4 (for inches) and t = time steps from 1 to
>>>>>>>>> 129.
>>>>>>>>>
>>>>>>>>> When I try to convert this to snow, I either get snow totals off
>>>>>>>>> the charts 100+ inches or I get rolling precip that moves and does not
>>>>>>>>> accumulate.  I am using a 10:1 ratio with csnowsfc to detect snow where 1
>>>>>>>>> is snow, 0 is no snow.
>>>>>>>>>
>>>>>>>>> I started with the equation above:
>>>>>>>>>
>>>>>>>>> if t=1|t=2
>>>>>>>>> 'display 'asnow'
>>>>>>>>> else
>>>>>>>>> 'define preciptot = sum(precip,t-1,t+0) - sum(precip,t-1,t-1)'
>>>>>>>>> 'display preciptot * 10'
>>>>>>>>> endif
>>>>>>>>>
>>>>>>>>> Where:
>>>>>>>>> asnow = apcpsfc * csnowsfc * 10.
>>>>>>>>>
>>>>>>>>> The above equation shows a 10:1 ratio for all of the precip.  Once
>>>>>>>>> I attempt to multiply by * csnowsfc (in the else section), I lose total
>>>>>>>>> precip and instead see moving precip over time.  I tried it in the sum
>>>>>>>>> equations multiplying precip * csnowsfc; then just in the display preciptot
>>>>>>>>> * csnowsfc * 10.  I even tried to use the same sum equation for precip
>>>>>>>>> substituting precip for csnowsfc, but none of these attempts worked to show
>>>>>>>>> total snow accumulations from t=1 to 129. I just can't get it to work with
>>>>>>>>> GFS.  All other models like nam, hrrr, use a nice simple sum equation:
>>>>>>>>> 'display sum(asnow,t=1,t='t')'.
>>>>>>>>>
>>>>>>>>> Any help here would be really appreciative.
>>>>>>>>>
>>>>>>>>> Sincerely,
>>>>>>>>>
>>>>>>>>> Jeff C
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> _______________________________________________
>>>>>>>>> gradsusr mailing list
>>>>>>>>> gradsusr at gradsusr.org
>>>>>>>>> http://gradsusr.org/mailman/listinfo/gradsusr
>>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> --
>>>>>>>> Jeff Duda, Research Scientist
>>>>>>>> University of Colorado Boulder
>>>>>>>> Cooperative Institute for Research in Environmental Sciences
>>>>>>>> NOAA/OAR/ESRL/Global Systems Laboratory
>>>>>>>> Boulder, CO
>>>>>>>> _______________________________________________
>>>>>>>> gradsusr mailing list
>>>>>>>> gradsusr at gradsusr.org
>>>>>>>> http://gradsusr.org/mailman/listinfo/gradsusr
>>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> gradsusr mailing list
>>>>>>> gradsusr at gradsusr.org
>>>>>>> http://gradsusr.org/mailman/listinfo/gradsusr
>>>>>>>
>>>>>>
>>>>>>
>>>>>> --
>>>>>> Jeff Duda, Research Scientist
>>>>>> University of Colorado Boulder
>>>>>> Cooperative Institute for Research in Environmental Sciences
>>>>>> NOAA/OAR/ESRL/Global Systems Laboratory
>>>>>> Boulder, CO
>>>>>> _______________________________________________
>>>>>> gradsusr mailing list
>>>>>> gradsusr at gradsusr.org
>>>>>> http://gradsusr.org/mailman/listinfo/gradsusr
>>>>>>
>>>>> _______________________________________________
>>>>> gradsusr mailing list
>>>>> gradsusr at gradsusr.org
>>>>> http://gradsusr.org/mailman/listinfo/gradsusr
>>>>>
>>>>
>>>>
>>>> --
>>>> Jeff Duda, Research Scientist
>>>> University of Colorado Boulder
>>>> Cooperative Institute for Research in Environmental Sciences
>>>> NOAA/OAR/ESRL/Global Systems Laboratory
>>>> Boulder, CO
>>>> _______________________________________________
>>>> gradsusr mailing list
>>>> gradsusr at gradsusr.org
>>>> http://gradsusr.org/mailman/listinfo/gradsusr
>>>>
>>> _______________________________________________
>>> gradsusr mailing list
>>> gradsusr at gradsusr.org
>>> http://gradsusr.org/mailman/listinfo/gradsusr
>>>
>>
>>
>> --
>> Jeff Duda, Research Scientist
>> University of Colorado Boulder
>> Cooperative Institute for Research in Environmental Sciences
>> NOAA/OAR/ESRL/Global Systems Laboratory
>> Boulder, CO
>> _______________________________________________
>> gradsusr mailing list
>> gradsusr at gradsusr.org
>> http://gradsusr.org/mailman/listinfo/gradsusr
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://gradsusr.org/pipermail/gradsusr/attachments/20220104/1a83702b/attachment-0001.html>

More information about the gradsusr mailing list