===============================================================
Earlier I had sent this mail by choosing "reply option" instead
of composing a new one...SORRY for the inconvenience caused.
===============================================================
Dear Ferreters,
Again there is a question regarding climatology. There are very
few mails in the archive regarding climatology and those are not helpfull
for my problem(?).
1. How to define the time axis to find the 'exact' climatology
for
a. daily data with one or more leap year in the data range
b. hourly data (say 6 hr) data with one or more leap year in
the data range
Regridding to non-leap time axis is an answer to above two
question. But in that case we will be forced to define a
new variable and use some kind of regridding utilities like
@ASN. In this case data for Feb-29 th of the leap year will
be assigned to Mar-01 st, which i don't like. The idea is that
data for feb-28th of all year should be averaged for getting
climatology for the day feb-28, and should be like that for
all days in the year (skipping feb-29 th of all leap years from
calculation)
If you are interested in this question please have a trial with
theattached jnl file.
2. Suppose we want to find the climatology explicitly. We can find
the average within a repeat loop but how to save this value
to an array (1d array of time) so that i can use it outside
the repeat loop ?
For example consider two year case and suppose we can skip
Feb-29th of every leap year from averaging;
repeat/k=1:365 ( ;\
;\
let mean = (var[t=`day`-`month`-`year_1`] +
var[t=`day`-`month`-`year_2`])/2 ;\
)
in this case, out side the repeat loop if we are listing 'mean',
we will get only one value for the DEC-31st. Definig a
climatological time axis and assigning that time-grid to the
variable 'mean'( before repeat loop) will not help. Ofcourse we
can write 'mean' to a file inside the repeat loop.
Does there is any other easy way ?
thanks in advance
Jaison Kurian
! climatology and leap year
! -------------------------
! let us find the daily climatology for some dummy data
! defining daily time axis for 1991-1992 (with leap year 1992) and creating
! dummy data
define axis/t="01-jan-1991":"31-dec-1992":1/unit=days tdly
let u_dly = COS(t[gt=tdly]/100) * SIN(t[gt=tdly]/100)+10
message
DEFINE AXIS/T=0:365.2425/EDGES/NPOINTS=365/T0="1-JAN-0001"/UNITS=DAYS/MODULO tdly_clim
let u_dly_clim = u_dly[gt=tdly_clim@MOD]
! Let us calculate the climatology by regridding to tdly_clim and
! compare it with explicitly calculated value
message
! Here it is seen that there is a shift in the climatology by one day, i.e.
! ( u_dly[t=02-jan-1991] + u_dly[t=02-jan-1992] )/2 =
! u_dly_clim[t=01-jan-0001]
! Let us see the Error in the climatology calculated by regridding.
! NOTE: average of 02-jan-1991 & 02-jan-1992 will be compared with
! 01-jan-0001 and so on
! WITH FEB-29 : data of Feb-29th-1991 will be averaged with march-01-1992
! WITHOUT FEB-29 : averaging is for same month same day's data only
! listing is as follows:
! "CLIM_DAY" ERROR ERROR
! (WITH FEB-29) (WITHOUT FEB-29)
message
\ cancel mode verify
repeat/k=1:364 ( ;\
let l_with_feb29 = `k+365` ;\
let l_no_feb29 = if `k+365 GE 425` then `k+366` else `k+365` ;\
let mean_with_f29 = ( u_dly[l=`k+1`] + u_dly[l=`l_with_feb29+1`] ) /2 ; \
let mean_no_f29 = ( u_dly[l=`k+1`] + u_dly[l=`l_no_feb29+1`] ) /2 ; \
let error_with_f29 = mean_with_f29 - u_dly_clim[l=`k`] ; \
let error_no_f29 = mean_no_f29 - u_dly_clim[l=`k`] ; \
if `k EQ 61 OR k EQ 277` then message ERROR > 0 ;\
list/nohead/format=(f9.4,3x,2f11.4)/quiet `k`, error_with_f29, error_no_f29 ;\
)
message
\ set mode verify
! without Feb-29 : Error is not zero from clim_day 60 to 275....!!!!
! with Feb-29 : Error is not zero from clim_day 275 onwards....!!!!
! what causes the shift by one day ?
message
! Let us use a different definition for climatological axis
DEFINE AXIS/T="1-JAN-0001":"31-DEC-0001":1/UNITS=days/MODULO tdly_clim
let u_dly_clim = u_dly[gt=tdly_clim@MOD]
! in this case there is no systematic shift as in previous case
! listing is as follows:
! CLIM_DAY ERROR ERROR
! (WITH FEB-29) (WITHOUT FEB-29)
message
\ cancel mode verify
repeat/k=1:364 ( ;\
let l_with_feb29 = `k+365` ;\
let l_no_feb29 = if `k+365 GE 425` then `k+366` else `k+365` ;\
let mean_with_f29 = ( u_dly[l=`k`] + u_dly[l=`l_with_feb29`] ) /2 ; \
let mean_no_f29 = ( u_dly[l=`k`] + u_dly[l=`l_no_feb29`] ) /2 ; \
let error_with_f29 = mean_with_f29 - u_dly_clim[l=`k`] ; \
let error_no_f29 = mean_no_f29 - u_dly_clim[l=`k`] ; \
list/nohead/format=(f9.4,3x,2f11.4)/quiet `k`, error_with_f29, error_no_f29 ;\
)
message
\ set mode verify
! in this case Error > 0 for all days and for both cases
message
! Let us do the same exercise for HOURLY data
define axis/t="01-jan-1991:00:00":"31-dec-1992:18:00":6/unit=hours t6hr
let u_6hr = COS(t[gt=t6hr]/100) * SIN(t[gt=t6hr]/100) + 10
message
! defining climatological axis
DEFINE AXIS/T=0:`365.2425*24`/EDGES/NPOINTS=`365*4`/T0="1-JAN-0001:00:00"/UNITS=hours/MODULO t6hr_clim
let u_6hr_clim = u_6hr[GT=t6hr_clim@MOD]
! no systematic shift is observed in this case
! listing (for the first 300 points only) is as follows:
! CLIM_DAY ERROR ERROR
! (WITH FEB-29) (WITHOUT FEB-29)
\ cancel mode verify
message
repeat/k=1:300 ( ;\
let l_with_feb29 = `k+1460` ;\
let l_no_feb29 = if `k+1460 GE 1697` then `k+1464` else `k+1460` ;\
let mean_with_f29 = ( u_6hr[l=`k`] + u_6hr[l=`l_with_feb29`] ) /2 ; \
let mean_no_f29 = ( u_6hr[l=`k`] + u_6hr[l=`l_no_feb29`] ) /2 ; \
let error_with_f29 = mean_with_f29 - u_6hr_clim[l=`k`] ; \
let error_no_f29 = mean_no_f29 - u_6hr_clim[l=`k`] ; \
list/nohead/format=(f9.4,3x,2f11.4)/quiet `k`, error_with_f29, error_no_f29 ;\
)
message
\ set mode verify
! in this case Error > 0 for all days and for both cases
! Let us use a different definition for climatological axis
DEFINE AXIS/T="1-JAN-0001:00:00":"31-DEC-0001:18:00":6/UNITS=hours/MODULO t6hr_clim
let u_6hr_clim = u_6hr[GT=t6hr_clim@MOD]
! no systematic shift is observed in this case
! listing (for the first 300 points only) is as follows:
! CLIM_DAY ERROR ERROR
! (WITH FEB-29) (WITHOUT FEB-29)
\ cancel mode verify
message
repeat/k=1:300 ( ;\
let l_with_feb29 = `k+1460` ;\
let l_no_feb29 = if `k+1460 GE 1697` then `k+1464` else `k+1460` ;\
let mean_with_f29 = ( u_6hr[l=`k`] + u_6hr[l=`l_with_feb29`] ) /2 ; \
let mean_no_f29 = ( u_6hr[l=`k`] + u_6hr[l=`l_no_feb29`] ) /2 ; \
let error_with_f29 = mean_with_f29 - u_6hr_clim[l=`k`] ; \
let error_no_f29 = mean_no_f29 - u_6hr_clim[l=`k`] ; \
list/nohead/format=(f9.4,3x,2f11.4)/quiet `k`, error_with_f29, error_no_f29 ;\
)
message
\ set mode verify
! in this case also Error > 0 for all days and for both cases
! So ... how to handle leap year for getting accurate climatology ? ????