;sample code to plot the annual average temperature time series
;for the 9 grid cells around the city of Riverside, CA.
pro climatology_501_timeseries_bias_correction
;where our downloaded .NC files are at
inputFolder = '/path/to/your/downloaded/data/'
;a list of the files we want to plot
files = inputFolder+['Riverside_NCEP.nc', $
'Riverside_ECH5.nc',$
'Riverside_GENMOM.nc',$
'Riverside_GFDL.nc']
; a list of model labels
modelLabels = ['NOAA NCEP',$
'MPI ECHAM5',$
'USGS GENMOM',$
'GFDL CM2.0']
;open the NCEP netcdf file and retrieve data
ncid0=ncdf_open(files[0], /NOWRITE)
NCDF_VARGET, ncid0, "time",ncep_time
;the temperature field is 3x3xN
NCDF_VARGET, ncid0, "TA",ncep_temperature
;close the netcdf
NCDF_CLOSE, ncid0
;average the 3x3 grid cells into one time series
ncep_temperature = mean(mean(ncep_temperature,DIMENSION=1),DIMENSION=1)
;convert model times to Julian day using 01-01-1900 as the base date
ncep_julian_date = ncep_time + JULDAY(1,1,1900)
;convert julian date into months, days, years
CALDAT, ncep_julian_date, ncep_months, ncep_days, ncep_years
;find the array indexes that corresponds to the months between 1980-1999
indexes = where(ncep_years GE 1980 AND ncep_years LE 1999)
;average the monthly temperature between 1980-1999 to give us a
;20-year climatology
ncep_annual_climatology_temperature = mean(ncep_temperature[indexes])
;setup graphics window size and title
w = window(DIMENSIONS=[800,800], $
WINDOW_TITLE="Riverside, CA Annual Temperature" )
;loop over each file to plot
for fileID=0,N_ELEMENTS(files)-1 do begin
;open the netcdf file and retrieve data
ncid0=ncdf_open(files[fileID], /NOWRITE)
NCDF_VARGET, ncid0, "time",time
;the temperature field is 3x3xN
NCDF_VARGET, ncid0, "TA",temperature
;close the netcdf
NCDF_CLOSE, ncid0
;average the 3x3 grid cells into one time series
temperature = mean(mean(temperature,DIMENSION=1),DIMENSION=1)
;keep a variable for how many time steps are in this file
numberOfMonths = N_ELEMENTS(time)
;reform data from N to [12, N/12] so we have each month in a matrix
temperature_by_month = reform(temperature, [ 12, numberOfMonths/12])
;average the first dimension of the matrix so w have an annual
;average temperature of length N/12
annual_average_temperature = mean(temperature_by_month, DIMENSION=1)
;similar to temperature, reform the time stamps to match the same
;[12, N/12] shape so that the time stamps match the original data
time_by_month = reform(time, [12, numberOfMonths/12])
;find the middle of the year. Should be around the beginning of July.
annual_average_time = mean(time_by_month, DIMENSION=1)
;convert model times to Julian day using 01-01-1900 as the base date
julian_date = annual_average_time + JULDAY(1,1,1900)
;convert julian date into months, days, years
CALDAT, julian_date, months, days, years
;find the array indexes that corresponds to the months between 1980-1999
indexes = where(years GE 1980 AND years LE 1999)
;average the monthly temperature between 1980-1999 to give us a
;20-year climatology
annual_climatology_temperature = mean(annual_average_temperature[indexes])
;calculate the difference (bias) between the current model climatology
;and the NCEP 1980-1999 climatology
biasValue = annual_climatology_temperature - ncep_annual_climatology_temperature
;display the bias in the IDL terminal
print, modelLabels[fileID], biasValue
;remove the bias from the annual average temperature
bias_corrected_temperature = annual_average_temperature - biasValue
;-----OPTIONAL CODE TO REMOVE GAP-----
delta_years = deriv(years)
gap_index = where(delta_years gt 1, count)
if(count gt 0 ) then begin
firstGap = gap_index[0]
;inject a NAN between the present data and the future data so that a
;straight line is not drawn by IDL
bias_corrected_temperature = [bias_corrected_temperature[0:firstGap], !VALUES.F_NAN, bias_corrected_temperature[firstGap+1:-1]]
julian_date = [julian_date[0:firstGap], !VALUES.F_NAN, julian_date[firstGap+1:-1]]
endif
;plot
p = plot(julian_date, bias_corrected_temperature, $
TITLE=modelLabels[fileID],LAYOUT=[2,2,fileid+1], /CURRENT, $
XRANGE=[JULDAY(1,1,1968),JULDAY(1,1,2100)], XTICKUNITS="Years", $
YRANGE=[12,22], YTITLE="Riverside, CA Temperature (!Eo!NC)",$
FONT_SIZE=12)
endfor
end