;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