;---------------------------------------------------------------------------- ; $Id: calcanom_spl.pro,v 1.2 2003/04/04 03:58:54 wblin Exp wblin $ ;+ ; Name: ; CALCANOM ; ; Purpose: ; Calculate anomalies given a continuous x,y,t hyperslab of data (i.e. ; timeseries of x-y slabs), by calculating monthly climatology and ; fitting the climatology with a spline fit. ; ; NB: This procedure assumes a year is 365 days long. Input time ; interval can be regular or irregular, although if input is CCM3 ; output monthly means, the keyword Ccm3_Moly_Means below needs to ; be set. Procedure has been tested with daily data, CCM3 monthly ; means, and 12 hourly data. ; ; Category: ; CCM3 Tools, Statistics. ; ; Calling Sequence: ; CALCANOM, in_time, in_date, io_data ; ; File I/O: ; None. ; ; Entry and Exit States: ; N/A. ; ; Input Parameters: ; in_time Vector of time at each timestep. Must be same size as ; the last dimension of io_data and monotonically ; increasing. Assumed to be in days. ; ; in_date Vector of integer values of the date at each in_time ; element in format yyyymmdd (with leading zeros absent). ; Must be same size as the last dimension of io_data and ; monotonically increasing in time. ; ; Input and Output Parameter: ; io_data Hyperslab (dimensioned x,y,time) of data. Overwritten ; with calculated anomalies and returned. Retains input ; type. Will work for array dimensioned [1,1,time]. ; ; Keywords (optional unless otherwise noted): ; BASE_DATE The date (in same format as in_date) that in_time=0 ; occurs at. Default value is 1st element of in_date. ; Integer. ; ; BASE_SEC The number of seconds on the date specified by Base_Date ; that in_time=0 occurs at. Default value is 0. Floating. ; ; CCM3_MOLY_MEANS If keyword set true, input data in_data consists of ; monthly means generated by CCM3. The unique quirk ; of this data is that in_date is written the day after ; the month the mean applies to (i.e. a Sep 00 mean has ; date 1001), which is compensated for by this function. ; ; Revision History: ; - 30 Aug 2001: Orig. ver. Johnny Lin, CIRES/University of Colorado. ; Email: air_jlin@yahoo.com. Passed moderately reasonable tests. ; - 19 Sep 2001: FOUND MAJOR BUG! ANY PRIOR VERSIONS OF THIS ROUTINE ; GIVES INCORRECT RESULTS. Basically, I found out in previous versions ; the spline-fit climatology is output, not the anomalies. Passed ; passably reasonable tests. ; - 6 May 2002: Fixed small bug that causes program to stop executing if ; in_time[0] is more than ~1095 days. Passed minimally passably ; reasonable tests. ; ; Notes: ; - Written for IDL 5.3. ; - This procedure cannot calculate anomalies for a timeseries of ; less than one full year. ; - Note that time and date generally proceed differently. For instance, ; if midnight beginning of day 1 is time=0, then midnight beginning of ; day 2 (i.e. midnight end of day 1) is time=1. Then time=2.5 would be ; noon of day 3. ; - No procedures called with _Extra keyword invoked. ; - User-written procedures called: CALCCLIM. ; - No common blocks are used in this program. ;- ; Copyright (c) 2001 Johnny Lin. For licensing and contact information ; see http://www.johnny-lin.com/lib.html. ;---------------------------------------------------------------------------- PRO CALCANOM, in_time, in_date, io_data $ , BASE_DATE = in_base_date $ , BASE_SEC = in_base_sec $ , CCM3_MOLY_MEANS = ccm3_moly_means $ , _EXTRA = extra ; -------------------- Error Check and Parameter Setting -------------------- ON_ERROR, 0 FORWARD_FUNCTION CALCCLIM, MEAN if (N_PARAMS() ne 3) then MESSAGE, 'error--bad number of params' time = in_time ;- protect input only parameters (rem that date = in_date ; time is assumed to be in days) sd_data = SIZE(io_data) ;- read and set array dims. if (sd_data[0] ne 3) then $ MESSAGE, 'error--data must be 3-D' NX = sd_data[1] NY = sd_data[2] NT = sd_data[3] if (N_ELEMENTS(in_base_date) ne 0) then $ ;- set base_date base_date = in_base_date $ else $ base_date = date[0] if (N_ELEMENTS(in_base_sec) ne 0) then $ ;- set base_sec base_sec = FLOAT(in_base_sec) $ else $ base_sec = 0.0 if ((SIZE(base_date, /Type) ne 2) and $ ;- check if is integer type (SIZE(base_date, /Type) ne 3)) then $ MESSAGE, 'error--not integer' if ((SIZE(date, /Type) ne 2) and $ (SIZE(date, /Type) ne 3)) then $ MESSAGE, 'error--not integer' if ((time[NT-1L]-time[0]+1.0) lt 365.0) then $ ;- check if timeseries is MESSAGE, 'error--less than 1 full yr data' ; more than 1 full year if (N_ELEMENTS(time) ne NT) then MESSAGE, 'error--bad dim' ;- check time dim if (N_ELEMENTS(date) ne NT) then MESSAGE, 'error--bad dim' ; consistent tmp = WHERE(date lt 101, count) ;- check no dates prior if (count gt 0) then MESSAGE, 'error--bad date' ; to 1 Jan 0000. ; ---------------------- Information About the Calendar --------------------- yr = date/10000L ;- vectors of the years, months, and days mo = (date-(yr*10000L))/100L ; of each elem. in date dy = date-(yr*10000L)-(mo*100L) dys_in_mo = [ 31, 28, 31, 30, 31, 30 $ ;- vector of a yr of days in each , 31, 31, 30, 31, 30, 31 ] ; month tdy_mid_mo = FLOAT(dys_in_mo) / 2.0 ;- time (in day) in each month ; that is the middle of the mo. base_yr = base_date/10000L ;- extract base yr, base_mo = (base_date-(base_yr*10000L))/100L ; mo., and day from base_dy = base_date-(base_yr*10000L)-(base_mo*100L) ; base_date ; ------------- Rough Check That Base Date/Sec Is Not Too Wrong ------------- ; ; Variable Descriptions: ; base_tdy_s0 Time of the base_date and base_sec in days since ; 1 Jan 0000 (i.e. 1-1-00), 00:00:00. ; date_s0 Time of date in days since 1-1-00, 00:00:00. ; diff_sec Diff. in sec between time_s0 and date_s0. ; tdy_of_yr_mo1 Time (in days since 1 Jan 0:00) that the 1st day of ; each month falls. ; time_s0 Time of time in days since 1-1-00, 00:00:00. ; Algorithm: For the each value of time and date, calculate what the ; values should be wrt to a fixed datum (1-1-00, 00:00:00). The diff. ; between the two should be less than 1 day, for each element. This ; routine operates on time and date as output by CCM3 (i.e. no adjustments ; for the "CCM3 monthly means" quirk). tdy_of_yr_mo1 = [ 0, 31, 59, 90, 120, 151 $ , 181, 212, 243, 273, 304, 334 ] base_tdy_s0 = (base_yr*365.0) + tdy_of_yr_mo1[base_mo-1] $ + (base_dy-1) + (base_sec/86400.0) time_s0 = time + base_tdy_s0 date_s0 = (yr*365.0) + tdy_of_yr_mo1[mo-1] + (dy-1.0) diff_sec = ROUND( (time_s0 - date_s0)*86400.0 ) tmp = WHERE(ABS(diff_sec) ge 86400, count) if (count gt 0) then begin MESSAGE, 'error--bad base keywords. ' $ + 'num. bad diff_sec = ' + STRTRIM(STRING(count),2) endif ; ------------ Adjust date and time For CCM3 Monthly Means Case ------------- ; ; Algorithm: time is moved back 1/2 month. I don't think I need to ; adjust date b/c CALCCLIM will automatically adjust for date as is, and ; the components of date are not used after this section. See description ; of Ccm3_Moly_Means keyword to see why time needs to be adjusted for this ; case. if (KEYWORD_SET(ccm3_moly_means) eq 1) then begin PRINT, 'CALCANOM: Note--Input data is CCM3 monthly means' prev_mo = mo-1 ;- previous month prev_mo_idx = prev_mo-1 ;- index of previous month dec_pts = WHERE(prev_mo_idx eq -1, count) ;- if previous month index =-1 if (count gt 0) then $ ; it means it is December; prev_mo_idx[dec_pts] = 11 ; chg. index to 11 time = time - tdy_mid_mo[prev_mo_idx] ;- subtract from time 1/2 of ; previous month endif ; ---------------------- Calculate Monthly Climatology ---------------------- clim_tmp = CALCCLIM( date, io_data $ ;- moly climatology , Ccm3_Moly_Means=ccm3_moly_means ) ;- reorder time dimens. so the values corresponding to month equal to ; base_mo is the 1st element in the time dimens. ("bms" means "base_mo ; start") clim_bms = SHIFT( clim_tmp, 0,0,-(base_mo-1L) ) ;+ reorder climatology tdy_mid_mo_bms = SHIFT( tdy_mid_mo, -(base_mo-1L) ) ;+ reorder tdy_mid_mo ;- calc. at what time the middle of the month (which is when we take the ; climatology as "being") occurs tdy_clim_bms = FLTARR(12) tdy_clim_bms[0] = tdy_mid_mo_bms[0] $ - (base_dy - 1.0) - (base_sec/86400.0) Ntcb = N_ELEMENTS(tdy_clim_bms) ;+ number elements in tdy_clim_bms for imo=1L,Ntcb-1L do begin tdy_clim_bms[imo] = tdy_clim_bms[imo-1] $ + tdy_mid_mo_bms[imo-1] $ + tdy_mid_mo_bms[imo] endfor ;- expand tdy_clim_bms from -(2-3) years from the beginning of vector time ; to an extra 3-4 years after the end of vector time to enable a spline ; fit to be calculated; expand_factor is used to count how many copies ; of a year to make to add to the beginning and end of time. the total ; number of such copies is factor_tot; tdy_clim is tdy_clim_bms expanded, ; and the corresponding climatology array to those times is clim. tdy_clim = tdy_clim_bms expand_factor = CEIL( (time[NT-1L]-tdy_clim_bms[Ntcb-1L]) $ ;+ expand after / 365.0 ) + 4 ; time vector if (expand_factor lt 0) then MESSAGE, 'error--strange sit.' for ifact=1L,expand_factor do begin tdy_clim = [ TEMPORARY(tdy_clim) $ , tdy_clim_bms + (ifact*365.0) ] endfor factor_tot = expand_factor ;> initialize tot. num factors expand_factor = CEIL( (tdy_clim_bms[0]-time[0]) $ ;+ expand before / 365.0 ) + 3 ; time vector if (expand_factor lt 3) then expand_factor=3 for ifact=1L,expand_factor do begin tdy_clim = [ tdy_clim_bms - (ifact*365.0) $ , TEMPORARY(tdy_clim) ] endfor factor_tot = factor_tot + expand_factor ;> expand tot. num factors clim = clim_bms for ifact=1L,factor_tot do begin clim = [ [[TEMPORARY(clim)]], [[clim_bms]] ] endfor sd_clim = SIZE(clim) if (sd_clim[3] ne N_ELEMENTS(tdy_clim)) then $ ;+ check num elements ok MESSAGE, 'error--bad expansions' ; ------ Spline Interpolate from Climatology to Timeseries and Subtract ----- ; ; Algorithm: Cycle through each spatial point. Select subrange that ; runs from the first annual climatological max. to the last annual ; climatological max. in the expanded timeseries. This subrange should ; be a superset in time of the vector time. Spline fit at the points ; defined by the vector time. Then subtract spline fit climatology to ; obtain the anomaly. PRINT, 'CALCANOM: Spline interpolation and calculate anomaly' for j=0L,NY-1L do begin for i=0L,NX-1L do begin clim_ij = REFORM(clim[i,j,*]) max_pts = WHERE(clim_ij eq MAX(clim_ij), count) if (count gt 0) then begin max1 = max_pts[0] ;- find loc. 1st global max. pt. max2 = max_pts[count-1L] ;- find loc. last global max. pt. endif if (tdy_clim[max1] ge time[0]) then MESSAGE, 'error--bad subrange' if (tdy_clim[max2] le time[NT-1L]) then MESSAGE, 'error--bad subrange' spl_type = SPL_INIT( tdy_clim[max1:max2] $ ;- spline w/ 1st , clim_ij[max1:max2] $ ; and last pt. , Yp0=0.0, Ypn_1=0.0 $ ; slope = 0 to , /Double ) ; obtain spline spl_clim = SPL_INTERP( tdy_clim[max1:max2] $ ; fit climatol. , clim_ij[max1:max2] $ , spl_type, time $ , /Double ) io_data[i,j,*] = io_data[i,j,*] - spl_clim[*] ;- calc. anom. endfor endfor ; --------------------------- Return From Procedure ------------------------- RETURN END ; ===== end of procedure ===== ; ========== end file ==========