Re: [ferret_users] seasonal mask

[Jaison Kurian <jaison@xxxxxxxxxxxxxxxxxx>]

Hi Prasad,
            If your data is for the period 1948-1997, it doesn't include
years 2000 & 2005. Hence you are left with only 4 "strong" years.

 The "easy" method will be extracting the 4 years data to 4 different
variables and then finding the mean. Like

    define axis/t=01-JAN-0001:31-DEC-0001:1/units=days/T0=31-DEC-0000/\
                    MODULO  tclim
    let wind_1950_1 = wind[t=01-JAN-1950:31-DEC-1950]
    let wind_1950   = wind_1950_1[gt=tclim@ASN]
    .................................................
    .................................................
    let wind_1970_1 = wind[t=01-JAN-1970:31-DEC-1970]
    let wind_1970   = wind_1970_1[gt=tclim@ASN]
    let wind_strg_mean = ( wind_1950 + wind_1965 + wind_1966 + wind_1970 )/4

But as the number of years to be included in the calculation increases,
there is a possibility for running out of memory. Hence the neat way to
do this is to apply a mask over time axis and then find the climatology
using @MOD regridding (as you tried to do). But here the problem is that
if the calendar for the variable is not NOLEAP/360_DAY, then you will not
get an exact match with the above method. Hence it will be better if you
can regrid your data to a noleap calendar time axis before using the
masking method described below.

Have a look at the example below. If you need more explanation regarding
any step, please let me know.

Hope this helps,

Regards,

Jaison

!----------Example---------------------------------------------------------

! let us create a dummy variable

     define axis/t=01-JAN-1948:31-DEC-1997:1/units=days/\
                         T0=31-DEC-1947/CALENDAR=NOLEAP time
     let wind = SIN(t[gt=time]/100)

! get the "L" indices of "strong" days in the "strong" years in a text file

     cancel reg/all  ! keep it here itself

     let l_start = `wind,r=lstart`
     let l_end   = `wind,r=lend`

     let yr_strong = TSEQUENCE({1950,1965,1966,1970}) ! should be in the
                                                      ! ascending order.
     sp rm -f dummy_L.dat
     REPEAT/RANGE=1:`yr_strong,r=lend`:1/NAME=yy ( ;\
        let cyr = yr_strong[l=`yy`] ;\
        list/nohead/file=dummy_L.dat/format=(3x,f10.1)/quiet/append \
                   L[gt=wind,t=01-JAN-`cyr`:31-DEC-`cyr`]  ;\
      )

! read in this "L" indices and map it back to original time axis using @XACT

     let strong_days = `{spawn:"cat dummy_L.dat | wc -l"}`

     define axis/t=1:`strong_days`:1 tfile
     define grid/t=tfile gfile

     FILE/grid=gfile/var=lindx dummy_L.dat
     define axis/t/from_var tdum=lindx   ! use d=1 or d=2 stuff if needed
     let strg = t[gt=tdum]*0+100

     define axis/t=`l_start`:`l_end`:1 tabs

     let strg_map =  strg[gt=tabs@XACT]

! create a mask

     let strg_mask = IF strg_map THEN 1
     let strg_wind = wind[gt=tabs@ASN] * strg_mask

     plot  wind[gt=tabs@ASN], strg_wind

! find the mean and do a cross check

     define axis/t=01-JAN-0001:31-DEC-0001:1/units=days/T0=31-DEC-0000/\
                    MODULO/CALENDAR=NOLEAP  tclim

     define axis/t=1:365:1/MODULO tabs_clim

     let strg_wind_mean = strg_wind[gt=tabs_clim@MOD]
     let strg_wind_clim = strg_wind_mean[gt=tclim@ASN] ! the final variable

     list/t=31-JAN-0001  strg_wind_clim

     let dummy = (wind[t=31-JAN-1950] + wind[t=31-JAN-1965] + \
                  wind[t=31-JAN-1966] + wind[t=31-JAN-1970] )/4

     list dummy

! remove intermediate files

     sp rm -f dummy_L.dat

!----------End of Example----------------------------------------------------

On Wed, 26 Apr 2006, Prasad Thoppil (Forn-Natl) wrote:

> I know this is rather a simple question. Suppose I have 50 years
> (1948-1997) of daily time-series of wind speed, how do I find average
> only for selected years. That is out of 50 years I want to find average
> for 6 years such as
>
> Let strong = {1950,1965,1966,1970,2000,2005}
> Then apply a mask if year_of_wind EQ strong THEN wind
>
> Though I tried my logic, I could not get it right.
>
> Once I have a masked time-series, it is possible to find 6 year mean
> (daily) by
> def axis/t=1-jan-0000:31-dec-0000:1/unit=days taxis
> let strong_avg = wind[gt=taxis@MOD]
>
> Thanks
>
> Prasad