;OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO pro makeresp, filename, nsrc, ndat, respdx, resp ;OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO ; read in responses close, 60 openr, 60, filename readf, 60, ncomment, nresp ; currently 246 in the G matrix pdmrow = fltarr(nsrc-1) ; this will be 26 pdmall = fltarr(nresp,nsrc) ; all responses (246x27) str = ' ' for i = 1, ncomment do readf, 60, str for i = 0, fix(nresp)-1 do begin readf, 60, pdmrow for j = 0, fix(nsrc)-2 do pdmall (i,j) = pdmrow (j) pdmall(i,nsrc-1) = 1. endfor pdmall (nresp-1,nsrc-1) = 0. for i = 0, fix(ndat)-1 do begin for j = 0, fix(nsrc)-1 do begin resp (i,j) = pdmall(respdx(i)-1,j) endfor endfor Close, 60 end ;OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO pro bayesinv, src, csrc, dat, cdat, pdm ;OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO ns = n_elements(src) nd = n_elements(dat) isrc = src idat = dat resid = dat - pdm # src for i = 0, fix(nd)-1 do resid[i] = resid[i]/cdat[i,i]^0.5 for i = 0, fix(nd)-1 do begin for j = 0, fix(ns)-1 do begin pdm[i,j] = pdm[i,j]*csrc[j,j]^0.5/cdat[i,i]^0.5 endfor endfor svdc, pdm, spdm, lpdm, rpdm picard = (rpdm) # resid picard[*] = picard[*]*spdm[*]/(spdm[*]*spdm[*]+1) src = transpose( lpdm) # picard src = transpose(csrc)^0.5 # src + isrc ; now data picard[*] = picard[*]*spdm[*] dat = picard ## transpose(rpdm) - resid for i = 0, fix(nd)-1 do dat[i] = idat[i] + cdat[i,i]^0.5*dat[i] ; now covariance ; fix nonbayesian to bayesian singular values spdm[*] = spdm[*]/(spdm[*]^2+1)^0.5 for i = 0, fix(ns)-1 do lpdm[*,i] = lpdm[*,i]*spdm[*] csrc = transpose(csrc)^0.5#(identity(ns)-transpose(lpdm)#lpdm)#csrc^0.5 return end ;OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO ;OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO pro transcom ;OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO ;OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO ; IDL wrapper for transcom 3 inversions ; Peter Rayner, Kevin Gurney, Rachel Law, October 2000 ; Program reads the number of models, model names, prior fluxes, and ; observational data. Then for each model makes a response function and ; carries out an inversion finally calculating a series of output ; diagnostics. ;OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO ; ********************************************************************** ; open the control file and read the control parameters ; ********************************************************************** close, 10 openr, 10, 'control.dat' readf, 10, nmodels model_names = strarr(nmodels) dum = ' ' str = ' ' FOR i = 0, nmodels-1 DO BEGIN readf, 10, format='(a1,a100)',dum,str model_names[i] = strtrim(str) ; print, model_names[i] ENDFOR print,'control file read. The number of models is: ',nModels readf, 10, str readf, 10, str readf, 10, allpre readf, 10, str readf, 10, noocn readf, 10, str readf, 10, nobio readf, 10, str readf, 10, noocnbio IF allpre EQ 1.0 THEN BEGIN title = '22 regions, all pre-subtraction tracers present' outinv = 'outmod.allpre.dat' outcov = 'poscov.allpre.dat' outplt = 'outmod.allpre.ps' outgrp = 'outgrp.allpre.dat' outdff = 'outdff.allpre.dat' outred = 'outred.allpre.dat' outsum = 'outsum.allpre.dat' ENDIF IF noocn EQ 1.0 THEN BEGIN title = '22 regions, no oceanic pre-subtraction' outinv = 'outmod.noocnpre.dat' outcov = 'poscov.noocnpre.dat' outplt = 'outmod.noocnpre.ps' outgrp = 'outgrp.noocnpre.dat' outdff = 'outdff.noocnpre.dat' outred = 'outred.noocnpre.dat' outsum = 'outsum.noocnpre.dat' ENDIF IF nobio EQ 1.0 THEN BEGIN title = '22 regions, no biospheric pre-subtraction' outinv = 'outmod.nobiopre.dat' outcov = 'poscov.nobiopre.dat' outplt = 'outmod.nobiopre.ps' outgrp = 'outgrp.nobiopre.dat' outdff = 'outdff.nobiopre.dat' outred = 'outred.nobiopre.dat' outsum = 'outsum.nobiopre.dat' ENDIF IF noocnbio EQ 1.0 THEN BEGIN title = '22 regions, no oceanic/bio pre-subtraction present' outinv = 'outmod.noocnbiopre.dat' outcov = 'poscov.noocnbiopre.dat' outplt = 'outmod.noocnbiopre.ps' outgrp = 'outgrp.noocnbiopre.dat' outdff = 'outdff.noocnbiopre.dat' outred = 'outred.noocnbiopre.dat' outsum = 'outsum.noocnbiopre.dat' ENDIF ; ********************************************************************** ; read in prior fluxes ; ********************************************************************** close, 20 openr, 20, 'prior.flux.dat' print, 'opened priors file' readf, 20, ncomment,nsrc ; # of comments, # of src values (27) ; (26 + 1 offset) ; declare the input arrays nunks = nsrc ; # of unknowns (= # of src values) srcnm = strarr(nunks) ; region names src = fltarr(nunks) ; prior fluxes csrc = fltarr(nunks,nunks) ; prior flux covariance str = ' ' FOR i = 1, ncomment do readf, 20, str FOR i = 0, fix(nsrc)-1 DO BEGIN readf, format='(a15,2f13.6)', 20, str,x,y src[i] = x srcnm[i] = str csrc[i,i] = y^2 ENDFOR Close, 20 print,'the prior fluxes have been read' ; declare the output arrays src2 = fltarr(nunks) ; posterior fluxes for each model csrc2 = fltarr(nunks,nunks) ; posterior flux covariance for each model msrc2 = fltarr(nunks) ; the average (over models) of the posterior fluxes mcsrc2 = fltarr(nunks,nunks) ; the average (over models) of the posterior covariances. cmsrc2 = fltarr(nunks,nunks) ; the covariance of the average model estimates for each ; region (how much each varies from the mean estimate) ; turn off the pre-subtracted fields as instructed by the control file IF noocn EQ 1.0 THEN src(3) = 0.0 IF nobio EQ 1.0 THEN src(2) = 0.0 IF noocnbio EQ 1 THEN BEGIN src(3) = 0.0 src(2) = 0.0 ENDIF ; ********************************************************************** ; open and read the data file ; ********************************************************************** close, 30 openr, 30, 'datafile.dat' print, 'opened datafile' readf, 30, ncomment,ndat ; declare arrays dat = fltarr(ndat) ; data dat2 = fltarr(ndat) ; inversion-predicted data cdat = fltarr(ndat,ndat) ; prior data covariance (from GV 2000) cdat2 = fltarr(ndat,ndat) ; data covariance respdx= intarr(ndat) ; response function index datnm = strarr(ndat) ; station names datlat= fltarr(ndat) ; latitudes ; read in the observed concentrations at the stations str = ' ' FOR i = 1, ncomment do readf, 30, str sumy = 0 ; RSD summation FOR i = 0, fix(ndat)-1 DO BEGIN readf, format='(a12,i6,2f11.3,2f9.2)',30,str,icdat,x,y,lat,lon cdat [i,i] = y^2 dat [i] = x respdx[i] = icdat datnm [i] = str datlat[i] = lat if (i lt ndat-1) then sumy = sumy + y ENDFOR print, 'mean RSD: ',sumy/(ndat-1) Close, 30 print,'a station dataset has been read. Number of stations: ',ndat ; ********************************************************************** ; declare arrays, call the inversion subroutines looping over models ; ********************************************************************** pdm = fltarr(ndat,nunks) ; partial derivative matrix (77x27) all_src2 = fltarr(nmodels,nunks) ; posterior fluxes for all models all_csrc2 = fltarr(nmodels,nunks,nunks) ; posterior flux covariances for all models all_dat2 = fltarr(nmodels,ndat) ; inversion-predicted data ; across models FOR i = 0, fix(nmodels)-1 DO BEGIN src2 = src csrc2 = csrc dat2 = dat cdat2 = cdat makeresp, model_names[i], nsrc, ndat, respdx, pdm ; The following makes a flask G matrix ; If i eq 0 then begin ; openw, 99, 'Gmat.small.dat' ; for j = 0, (ndat-1) do printf, format='(27f8.4)', 99, pdm[j,*] ; close, 99 ; Endif bayesinv, src2, csrc2, dat2, cdat2, pdm all_src2[i,*] = src2 all_csrc2[i,*,*] = csrc2 all_dat2[i,*] = dat2 ENDFOR ; ********************************************************************** ; open and write the model/region output file ; ********************************************************************** close, 15 openw, 15, outinv form1 = string(format='("(i3,",i2,"f8.4)")',nmodels) form2 = string(format='("(i3,",i2,"f8.2)")',nmodels) printf, 15, 'The posterior sources for each model and basis function region' printf, 15, format='(i5,a1,i5)', nunks, ',', nmodels for i = 0, nsrc-2 do printf, format=form1, 15, i+1, all_src2[*,i] printf, format=form2, 15, nsrc, all_src2[*,nsrc-1] printf, 15 printf, 15, 'Posterior source standard deviation for each model and region' printf, 15, format='(i5,a1,i5)', nunks, ',', nmodels for i = 0, nsrc-1 do printf, format=form1, 15, i+1, all_csrc2[*,i,i]^0.5 ; now calculate averages of both posterior sources and covariances for i = 0, nunks-1 do begin msrc2[i] = mean (all_src2[*,i]) for j = 0, nunks-1 do begin mcsrc2[i,j] = mean (all_csrc2[*,i,j]) endfor endfor ; standard deviation and covariance of posterior sources. for i = 0, nunks-1 do begin for j = 0, nunks-1 do begin cmsrc2[i,j] = mean ((all_src2[*,i]-msrc2[i])*(all_src2[*,j]-msrc2[j])) if i eq j then begin if cmsrc2[i,j] lt 0 then begin print, 'Warning cmsrc2 < 0 for i,j: ',i,j cmsrc2[i,j] = 0 endif endif endfor endfor printf, 15 printf, 15, 'Mean model source, mean model source sd, sd of model means' printf, 15, format='(i5,a1,i5)', nunks, ',', 3 for i = 0, nsrc-2 do printf, 15, i+1, msrc2[i], mcsrc2[i,i]^0.5, $ cmsrc2[i,i]^0.5, format='(i3,3f9.4)' printf, 15, nsrc, msrc2[nsrc-1],mcsrc2[nsrc-1,nsrc-1]^0.5, $ cmsrc2[nsrc-1,nsrc-1]^0.5,format='(i3,f9.2,2f9.4)' Close, 15 ; write out the posterior source covariance matrix close, 25 openw, 25, outcov FOR i = 0, nmodels-1 DO BEGIN bigG = STRPOS(model_names[i],'G',/REVERSE_SEARCH) printf, 25, STRMID(model_names[i],10,bigG-11) FOR j = 0, nsrc-1 do printf, 25, all_csrc2[i,*,j], format='(27f8.4)' printf, 25 ENDFOR Close, 25 print, 'The model/region output files have been written' ; ********************************************************************** ; open and write the model/group output file ; ********************************************************************** readf, 10, ngroups close, 35 openw, 35, outgrp groupname = strarr(ngroups) g_msrc2 = fltarr(ngroups) ; model mean posterior source g_msrc2[*] = 0.0 ; for each regional group g_mcsrc2 = fltarr(ngroups) ; model mean posterior source g_mcsrc2[*] = 0.0 ; covariance for each regional group g_cmsrc2 = fltarr(ngroups) ; covariance of the model mean g_cmsrc2[*] = 0.0 ; posterior source for each regional group g_allsrc2 = fltarr(ngroups,nmodels); model posterior source for g_allsrc2[*,*] = 0. ; each regional group g_allcsrc2 = fltarr(ngroups,nmodels); model posterior source g_allcsrc2[*,*] = 0. ; cov for each regional group FOR n = 0, ngroups-1 DO BEGIN readf, 10, format='(a1,a13)',dum,str groupname[n] = strtrim(str) grouplen = 0 readf, 10, grouplen ; number of regions in the group grouplist = intarr(grouplen) ; array of the region numbers in group readf, 10, grouplist for i = 0, grouplen-1 do begin g_msrc2[n] = g_msrc2[n] + msrc2[grouplist[i]] g_allsrc2[n,*] = g_allsrc2[n,*] + all_src2[*,grouplist[i]] for j = 0, grouplen-1 do begin g_mcsrc2[n] = g_mcsrc2[n] + mcsrc2[grouplist[i], grouplist[j]] g_allcsrc2[n,*] = g_allcsrc2[n,*] + all_csrc2[*,grouplist[i], grouplist[j]] g_cmsrc2[n] = g_cmsrc2[n] + cmsrc2[grouplist[i],grouplist[j]] endfor endfor endfor form3 = string(format='("(",i2,"f7.2,a14)")',nmodels) form4 = string(format='("(",i2,"f7.2,a14)")',nmodels) printf, 35, 'Model posterior source for each regional group' printf, format='(i5,a1,i5)', 35, ngroups,',',nmodels For n = 0, ngroups-1 do printf, format=form3, 35, g_allsrc2[n,*], $ groupname[n] printf, 35 printf, 35, 'Model posterior source standard deviation for each regional group' printf, format='(i5,a1,i5)', 35, ngroups,',',nmodels For n = 0, ngroups-1 do printf, format=form4, 35, $ (g_allcsrc2[n,*])^0.5, $ groupname[n] printf, 35 printf, 35, 'Model mean source, model mean source sd, and sd of model mean' printf, format='(i5,a1,i5)', 35, ngroups,',',3 For n = 0, ngroups-1 do printf, format='(3f7.2,a14)', 35, g_msrc2[n],$ (g_mcsrc2[n])^0.5, $ (g_cmsrc2[n])^0.5,groupname[n] Close, 10 Close, 35 ; ********************************************************************** ; open and write the predicted versus measured and uncertainty reduction ; files ; ********************************************************************** close, 45 openw, 45, outdff printf, 45, 'Data mismatch for each model and site (1st col: latitude)' printf, format='(i5,a1,i5)', 45, ndat,',',nmodels sumdat = fltarr(nmodels) ; sum of (mismatch^2 divided by the prior data cov.) sumdat = 0 form5 = string(format='("(",i2,"f8.3,a13)")',nmodels+1) for i = 0, ndat-1 do begin if (cdat[i,i] lt 99) then begin printf, format=form5, 45, datlat[i],all_dat2[*,i]-dat[i], $ datnm[i] sumdat = sumdat + (all_dat2[*,i]-dat[i])^2.0/cdat[i,i] endif endfor Close, 45 close, 55 openw, 55, outred form7 = string(format='("(i4,",i2,"f8.1)")',nmodels) printf, 55, 'Reduction in uncertainty for land/ocean regions' printf, format='(i5,a1,i5)', 55, nunks-5,',',nmodels for i = 4, 25 do begin printf, format=form7, 55, i-3, (1.0-((all_csrc2[*,i,i]^0.5)/ $ (csrc[i,i]^0.5)))*100.0 endfor Close, 55 ; ********************************************************************** ; open and write the summary output file ; ********************************************************************** close, 65 openw, 65, outsum sumdat=sumdat/ndat ; average the sumdat ratio across stations ; 1 = mismatch distribution consistent with prior uncert. ; 0 = exact match to data. sumsrc = fltarr(nmodels) ; sum of((post minus prior flux)^2/prior source cov) sumsrc = 0 for i = 4, 25 do sumsrc = sumsrc+(all_src2[*,i]-src[i])^2.0/csrc[i,i] print, mean(sumdat), mean(sumsrc/ndat) sumsrc = sumdat + sumsrc/ndat ; a chi-squared measure tracname = ['fos90','fos95','Neutral NEP','Takahashi Ocean', $ 'L01: Bor N America','L02: Temp N America', $ 'L03: Trop America','L04: South America', $ 'L05: Northern Africa','L06: Southern Africa', $ 'L07: Boreal Asia','L08: Temperate Asia', $ 'L09: Tropical Asia','L10:Australasia','L11: Europe', $ 'O01: N Pacific','O02: W Pacific','O03: E Pacific', $ 'O04: S Pacific','O05: Northern Ocean','O06: N Atlantic',$ 'O07: Tropical Atlantic','O08: S Atlantic', $ 'O09: Southern Ocean','O10: Trop Indian Ocean', $ 'O11: S Indian Ocean'] takflux = fltarr(nunks) takflux[0:14] = -1*msrc2[0:14] takflux[15:25] = [-0.5016, 0.1523,0.4671,-0.2314,-0.4389,-0.2903, $ 0.1286,-0.1289,-0.8858,0.1182,-0.5556] printf, 65, 'Model average statistics:' printf, 65 printf, 65, ' prior posterior posterior posterior' printf, 65, ' tracer flux stdev flux stdev sd model means total ocean flux' printf, 65 for i = 0, nsrc-2 do printf, 65, tracname[i], src[i], csrc[i,i]^0.5, $ msrc2[i], mcsrc2[i,i]^0.5, $ cmsrc2[i,i]^0.5,msrc2[i]+takflux[i],$ format='(a23,2x,2f8.4,1x,2f8.4,4x,f7.4,9x,f7.4)' printf, 65, 'offset', src[nsrc-1], csrc[nsrc-1,nsrc-1]^0.5, $ msrc2[nsrc-1], mcsrc2[nsrc-1,nsrc-1]^0.5, $ cmsrc2[nsrc-1,nsrc-1]^0.5, $ format='(a23,2x,2f8.2,1x,f8.2,f8.4,3x,f8.4)' printf, 65 printf, 65 printf, 65, 'Model average statistics for regional lumping' printf, 65 printf, 65, ' Posterior' printf, 65, 'Regional group source stdev stdev of model means' printf, 65 For n = 0, ngroups-1 do printf, 65, groupname[n], g_msrc2[n], $ (g_mcsrc2[n])^0.5, $ (g_cmsrc2[n])^0.5, $ format='(a13,3x,2f8.4,8x,f8.4)' printf, 65 printf, 65 printf, 65, 'station/model average total data mismatch: ', $ mean(sumdat), format='(a43,f5.3)' printf, 65 printf, 65, 'average total data mismatch for each model' form6 = string(format='("(",i2,"f8.2)")',nmodels) printf, format=form6, 65, sumdat printf, 65 printf, 65 printf, 65, 'model average chi squared: ', mean(sumsrc), $ format='(a27,f5.3)' printf, 65 printf, 65, 'chi squared for each model' printf, format=form6, 65, sumsrc sum1 = 0 for i = 4, 25 do sum1 = sum1 + mcsrc2[i,i] + cmsrc2[i,i] printf, 65 printf, 65 printf, 65, 'Total unc. (avg model cov + cov of model means): ',sum1,$ format='(a49,f8.3)' Close, 65 ;==================================================================== ; Make an array of box-whisker plots of the results for each region ; write this array of 24 plots into a postscript file" SET_PLOT,'PS' !P.FONT = 0 DEVICE, filename=outplt, /helvetica, /color, /landscape ; Set a bunch of colors LOADCT, 13 TVLCT, r, g, b, /get r(255) = 255B g(255) = 255B b(255) = 255B TVLCT, r, g, b black = 0B white = 255B red = 254B blue = 100B nRegions = 24 !P.MULTI = [0,6,4] map = [4,19,14,2,10,15,5,20,8,3,11,16,6,21,9,24,12,17,7,22,23,25,13,18] reg = ['fos90','fos95','Neutral NEP','Takahashi Ocean', $ 'Boreal N America', 'Temperate N America','Tropical America', $ 'S America','N Africa','S Africa','Boreal Asia', $ 'Temperate Asia','Trop Asia','Australasia', $ 'Europe','N Pacific','W Pacific','E Pacific','S Pacific', $ 'Northern Ocean','N Atlantic','Tropical Atlantic','S Atlantic',$ 'Southern Ocean','Tropical Indian Ocean','S Indian Ocean'] DEVICE, FONT_SIZE=14 titlesize=1.5 textSize=0.65 errThick = 4 FOR k=0,nRegions-1 DO BEGIN ; loop over basis regions ; Draw and label the "frame" of the plot, with no data plotMax = max(all_src2(*,0:nunks-2)) + $ max(all_csrc2(*,0:nunks-2,0:nunks-2))^0.5 PLOT, [0.,nModels], [-plotMax,plotMax],TITLE=reg[map[k]], $ XTICKS=nModels+1, XTICKV=[0.,FINDGEN(nModels)+0.5,nModels], $ /NODATA, /XSTYLE, /ystyle, COLOR=black, $ XTICKNAME=[' ',STRCOMPRESS(STRING(INDGEN(nModels)+1)),' '] ; For each model result, draw a colored rectangle corresponding ; to the source strength, with a black error bar FOR j=0,nModels-1 DO BEGIN ; j is the model index number strength = all_src2[j,map[k]] ; source strength for model j boxX = [j,j,j+1,j+1] ; x coordinates of filled rectangles boxY = [strength, 0, 0, strength] ; y coordinates ; Use blue for sinks, red for sources IF (strength LT 0) THEN barColor=blue ELSE barColor=red ; Draw the colored rectangle POLYFILL, boxX, boxY, color=barColor ; Draw an outline of the rectangle in black OPLOT, boxX, boxY, color=black ; Draw an "error bar" with covariance above and below rectangle xErr = [j+0.5, j+0.5] err = all_csrc2[j,map[k],map[k]]^0.5 yErr = [strength+err, strength-err] oplot, xErr, yErr, color=1, thick=errThick ENDFOR ; next model (j) ENDFOR ; next region (k) ; Add explanatory text !P.MULTI=0 XYOUTS, .5, 1.03, title, ALIGNMENT=0.5, size=titleSize, /NORMAL XYOUTS, .5, -0.02, ALIGNMENT=0.5, size=textsize, /NORMAL, $ 'Annual Mean Inversion by Region and Model' XYOUTS, .5, -0.05, ALIGNMENT=0.5, size=textSize, /NORMAL, $ '1=mod.1 2=mod.2 3=mod.3 4=mod.4 5=mod.5 6=mod.6 7=mod.7 8=mod.8 9=mod.9 10=mod.10 11=mod.11 12=mod.12 13=mod.13 14=mod.14 15=mod.15' DEVICE, /CLOSE SET_PLOT, 'X' ; ====================================================================== end