[gradsusr] Total Snow Accumulation with GFS .25 Degree

[Jeff Chabot]

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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://gradsusr.org/pipermail/gradsusr/attachments/20220103/77eb25a3/attachment-0001.html>