pro average_fcor ; Purpose ; ------- ; Reads in monthly fcor output from the CSIRO AGCM, derives a climatology and ; writes the values to a netCDF file. ; ; History ; -------- ; 2002 Nov 20 Steve Phipps Original version ; 2003 May 29 Steve Phipps Modified to calculate DJF and JJA averages, and ; to read in correct latitude and longitude ; values, rather than derive approximate ones. ; Also modified for the 32-bit integer version of ; the model, although the program will still work ; OK with output from the previous version. ; 2003 Jun 23 Steve Phipps Modified for the addition of a new field ; (ATHFA) to the model output ; 2004 Sep 13 Steve Phipps Modified for the addition of seven new fields ; to the model output ; 2004 Sep 30 Steve Phipps Modified to read the latitudes and longitudes ; from csiro_ogcm_ts_landsea.nc, rather than ; csiro_landsea.nc; also modified for five-digit ; year numbers ; 2004 Nov 26 Steve Phipps Modified to calculate standard errors; DJF and ; JJA averages dropped. ; Define constants lon = 64 lat2 = 56 nmonth = 12 nvar = 14 days = [31.0, 28.0, 31.0, 30.0, 31.0, 30.0, 31.0, 31.0, 30.0, 31.0, 30.0, 31.0] ; Initialise variables sfluxes_ann = fltarr(lon, lat2, 2, nvar) sfluxes_clim = fltarr(lon, lat2, nmonth, 2, nvar) ; Get run details from user run = "" year1 = 0 year2 = 0 read, run, prompt="Run name [e.g. b44] " read, year1, prompt="Initial year for climatology [e.g. 6] " read, year2, prompt="Final year for climatology [e.g. 35] " ; Declare arrays to hold data nyear = year2 - year1 + 1 sfluxes = fltarr(lon, lat2, nmonth, nyear, nvar) sfluxes_in = dblarr(lon, lat2) nsteps = 0l ; Loop over years for year = year1, year2 do begin years = strtrim(string(year), 2) if (year lt 10000) then years = '0' + years if (year lt 1000) then years = '0' + years if (year lt 100) then years = '0' + years if (year lt 10) then years = '0' + years ; Loop over months for month = 1, 12 do begin months = strtrim(string(month), 2) if (month lt 10) then months = '0' + months ; Open file for unformatted read access filename = 'fcor' + years + months + '.' + run openr, 1, filename, /f77_unformatted ; Read data readu, 1, nsteps for i = 0, nvar-1 do begin readu, 1, sfluxes_in sfluxes(*, *, month-1, year-year1, i) = sfluxes_in endfor ; Close file close, 1 endfor endfor ; Derive climatologies ann_temp = fltarr(nyear) for k = 0, nvar-1 do begin for j = 0, lat2-1 do begin for i = 0, lon-1 do begin ; Monthly climatology for month = 0, 11 do begin temp = sfluxes(i, j, month, *, k) stats = moment(temp) sfluxes_clim(i, j, month, 0, k) = stats(0) sfluxes_clim(i, j, month, 1, k) = sqrt(stats(1) / float(nyear)) endfor ; Annual means ann_temp(*) = 0.0 for month = 0, 11 do begin ann_temp = ann_temp + days(month) * sfluxes(i, j, month, *, k) / 365.0 endfor stats = moment(ann_temp) sfluxes_ann(i, j, 0, k) = stats(0) sfluxes_ann(i, j, 1, k) = sqrt(stats(1) / float(nyear)) endfor endfor endfor ; Read latitude and longitude values from csiro_ogcm_ts_landsea.nc ncid = ncdf_open("csiro_ogcm_ts_landsea.nc") lonid = ncdf_varid(ncid, "lonts") latid = ncdf_varid(ncid, "latts") ncdf_varget, ncid, lonid, lons ncdf_varget, ncid, latid, lats ncdf_close, ncid ; Derive name of output file year1s = strtrim(string(year1), 1) year2s = strtrim(string(year2), 1) if (year1 lt 10000) then year1s = '0' + year1s if (year1 lt 1000) then year1s = '0' + year1s if (year1 lt 100) then year1s = '0' + year1s if (year1 lt 10) then year1s = '0' + year1s if (year2 lt 10000) then year2s = '0' + year2s if (year2 lt 1000) then year2s = '0' + year2s if (year2 lt 100) then year2s = '0' + year2s if (year2 lt 10) then year2s = '0' + year2s filename = "fcor_" + run + "_climats_" + year1s + "-" + year2s + ".nc" ; Save data mths = indgen(12) + 1 vid = intarr(nvar) vseid = intarr(nvar) vannid = intarr(nvar) vannseid = intarr(nvar) ncid = ncdf_create(filename) londid = ncdf_dimdef(ncid, "longitude", lon) latdid = ncdf_dimdef(ncid, "latitude", lat2) mthdid = ncdf_dimdef(ncid, "month", nmonth) lonvid = ncdf_vardef(ncid, "longitude", londid, /float) latvid = ncdf_vardef(ncid, "latitude", latdid, /float) mthvid = ncdf_vardef(ncid, "month", mthdid, /short) vid(0) = ncdf_vardef(ncid, "taux", [londid, latdid, mthdid], /float) vid(1) = ncdf_vardef(ncid, "tauy", [londid, latdid, mthdid], /float) vid(2) = ncdf_vardef(ncid, "athf", [londid, latdid, mthdid], /float) vid(3) = ncdf_vardef(ncid, "atsf", [londid, latdid, mthdid], /float) vid(4) = ncdf_vardef(ncid, "atsi", [londid, latdid, mthdid], /float) vid(5) = ncdf_vardef(ncid, "u3", [londid, latdid, mthdid], /float) vid(6) = ncdf_vardef(ncid, "athfa", [londid, latdid, mthdid], /float) vid(7) = ncdf_vardef(ncid, "dsis", [londid, latdid, mthdid], /float) vid(8) = ncdf_vardef(ncid, "dsisb", [londid, latdid, mthdid], /float) vid(9) = ncdf_vardef(ncid, "dsfw", [londid, latdid, mthdid], /float) vid(10) = ncdf_vardef(ncid, "atsf1", [londid, latdid, mthdid], /float) vid(11) = ncdf_vardef(ncid, "atsf2", [londid, latdid, mthdid], /float) vid(12) = ncdf_vardef(ncid, "atsf3", [londid, latdid, mthdid], /float) vid(13) = ncdf_vardef(ncid, "atsfold", [londid, latdid, mthdid], /float) vseid(0) = ncdf_vardef(ncid, "tauxse", [londid, latdid, mthdid], /float) vseid(1) = ncdf_vardef(ncid, "tauyse", [londid, latdid, mthdid], /float) vseid(2) = ncdf_vardef(ncid, "athfse", [londid, latdid, mthdid], /float) vseid(3) = ncdf_vardef(ncid, "atsfse", [londid, latdid, mthdid], /float) vseid(4) = ncdf_vardef(ncid, "atsise", [londid, latdid, mthdid], /float) vseid(5) = ncdf_vardef(ncid, "u3se", [londid, latdid, mthdid], /float) vseid(6) = ncdf_vardef(ncid, "athfase", [londid, latdid, mthdid], /float) vseid(7) = ncdf_vardef(ncid, "dsisse", [londid, latdid, mthdid], /float) vseid(8) = ncdf_vardef(ncid, "dsisbse", [londid, latdid, mthdid], /float) vseid(9) = ncdf_vardef(ncid, "dsfwse", [londid, latdid, mthdid], /float) vseid(10) = ncdf_vardef(ncid, "atsf1se", [londid, latdid, mthdid], /float) vseid(11) = ncdf_vardef(ncid, "atsf2se", [londid, latdid, mthdid], /float) vseid(12) = ncdf_vardef(ncid, "atsf3se", [londid, latdid, mthdid], /float) vseid(13) = ncdf_vardef(ncid, "atsfoldse", [londid, latdid, mthdid], /float) vannid(0) = ncdf_vardef(ncid, "tauxann", [londid, latdid], /float) vannid(1) = ncdf_vardef(ncid, "tauyann", [londid, latdid], /float) vannid(2) = ncdf_vardef(ncid, "athfann", [londid, latdid], /float) vannid(3) = ncdf_vardef(ncid, "atsfann", [londid, latdid], /float) vannid(4) = ncdf_vardef(ncid, "atsiann", [londid, latdid], /float) vannid(5) = ncdf_vardef(ncid, "u3ann", [londid, latdid], /float) vannid(6) = ncdf_vardef(ncid, "athfaann", [londid, latdid], /float) vannid(7) = ncdf_vardef(ncid, "dsisann", [londid, latdid], /float) vannid(8) = ncdf_vardef(ncid, "dsisbann", [londid, latdid], /float) vannid(9) = ncdf_vardef(ncid, "dsfwann", [londid, latdid], /float) vannid(10) = ncdf_vardef(ncid, "atsf1ann", [londid, latdid], /float) vannid(11) = ncdf_vardef(ncid, "atsf2ann", [londid, latdid], /float) vannid(12) = ncdf_vardef(ncid, "atsf3ann", [londid, latdid], /float) vannid(13) = ncdf_vardef(ncid, "atsfoldann", [londid, latdid], /float) vannseid(0) = ncdf_vardef(ncid, "tauxannse", [londid, latdid], /float) vannseid(1) = ncdf_vardef(ncid, "tauyannse", [londid, latdid], /float) vannseid(2) = ncdf_vardef(ncid, "athfannse", [londid, latdid], /float) vannseid(3) = ncdf_vardef(ncid, "atsfannse", [londid, latdid], /float) vannseid(4) = ncdf_vardef(ncid, "atsiannse", [londid, latdid], /float) vannseid(5) = ncdf_vardef(ncid, "u3annse", [londid, latdid], /float) vannseid(6) = ncdf_vardef(ncid, "athfaannse", [londid, latdid], /float) vannseid(7) = ncdf_vardef(ncid, "dsisannse", [londid, latdid], /float) vannseid(8) = ncdf_vardef(ncid, "dsisbannse", [londid, latdid], /float) vannseid(9) = ncdf_vardef(ncid, "dsfwannse", [londid, latdid], /float) vannseid(10) = ncdf_vardef(ncid, "atsf1annse", [londid, latdid], /float) vannseid(11) = ncdf_vardef(ncid, "atsf2annse", [londid, latdid], /float) vannseid(12) = ncdf_vardef(ncid, "atsf3annse", [londid, latdid], /float) vannseid(13) = ncdf_vardef(ncid, "atsfoldannse", [londid, latdid], /float) ncdf_attput, ncid, lonvid, "units", "degrees_east" ncdf_attput, ncid, latvid, "units", "degrees_north" ncdf_attput, ncid, mthvid, "units", "months" ncdf_control, ncid, /endef ncdf_varput, ncid, lonvid, lons ncdf_varput, ncid, latvid, lats ncdf_varput, ncid, mthvid, mths for i = 0, nvar-1 do begin ncdf_varput, ncid, vid(i), sfluxes_clim(*, *, *, 0, i) ncdf_varput, ncid, vseid(i), sfluxes_clim(*, *, *, 1, i) ncdf_varput, ncid, vannid(i), sfluxes_ann(*, *, 0, i) ncdf_varput, ncid, vannseid(i), sfluxes_ann(*, *, 1, i) endfor ncdf_close, ncid end