1 | pro regrid |
---|
2 | filename='/data/obolmd/CMIP6/CO2_EMISSIONS_48x36/flux_2D_CO2_1867.nc' |
---|
3 | print, filename |
---|
4 | NETCDFREAD,filename,'flux',flux,dimflux |
---|
5 | NETCDFREAD,filename,'lat',lat,dimlat0 |
---|
6 | NETCDFREAD,filename,'lon',lon,dimlon0 |
---|
7 | NETCDFREAD,filename,'TIME',time,dimtime0 |
---|
8 | dimlat=dimlat0(0) |
---|
9 | dimlon=dimlon0(0) |
---|
10 | dimtime=dimtime0(0) |
---|
11 | print, 'dim flux=', dimflux |
---|
12 | A = float([ [3./4., 1./8., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1./8.],$ |
---|
13 | [1./8., 3./4., 1./8., 0., 0., 0., 0., 0., 0., 0., 0., 0.],$ |
---|
14 | [0., 1./8., 3./4., 1./8., 0., 0., 0., 0., 0., 0., 0., 0.],$ |
---|
15 | [0., 0., 1./8., 3./4., 1./8., 0., 0., 0., 0., 0., 0., 0.],$ |
---|
16 | [0., 0., 0., 1./8., 3./4., 1./8., 0., 0., 0., 0., 0., 0.],$ |
---|
17 | [0., 0., 0., 0., 1./8., 3./4., 1./8., 0., 0., 0., 0., 0.],$ |
---|
18 | [0., 0., 0., 0., 0., 1./8., 3./4., 1./8., 0., 0., 0., 0.],$ |
---|
19 | [0., 0., 0., 0., 0., 0., 1./8., 3./4., 1./8., 0., 0., 0.],$ |
---|
20 | [0., 0., 0., 0., 0., 0., 0., 1./8., 3./4., 1./8., 0., 0.],$ |
---|
21 | [0., 0., 0., 0., 0., 0., 0., 0., 1./8., 3./4., 1./8., 0.],$ |
---|
22 | [0., 0., 0., 0., 0., 0., 0., 0., 0., 1./8., 3./4., 1./8.],$ |
---|
23 | [1./8., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1./8., 3./4.] ]) |
---|
24 | A = float([ [1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],$ |
---|
25 | [0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],$ |
---|
26 | [0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],$ |
---|
27 | [0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],$ |
---|
28 | [0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],$ |
---|
29 | [0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],$ |
---|
30 | [0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],$ |
---|
31 | [0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],$ |
---|
32 | [0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],$ |
---|
33 | [0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],$ |
---|
34 | [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],$ |
---|
35 | [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.] ]) |
---|
36 | fluxinit=flux |
---|
37 | fluxinit=flux |
---|
38 | flux_check=flux |
---|
39 | for lo=0,dimlon-1 do begin |
---|
40 | for la=0,dimlat-1 do begin |
---|
41 | flux_check(lo,la,*) = invert(A) ## transpose(fluxinit(lo,la,*)) |
---|
42 | endfor |
---|
43 | endfor |
---|
44 | m_bloq = make_array(dimlon,dimlat,12,value=0) ; Matrice booléenne "mois à bloquer ou non" |
---|
45 | if total(where(flux_check lt 0)) ne -1 then m_bloq(where(flux_check lt 0)) = 1 |
---|
46 | ; Correction/adaptation de la matrice S&Z en fonction du masque booléen m_bloq |
---|
47 | for lo=0,dimlon-1 do begin |
---|
48 | for la=0,dimlat-1 do begin |
---|
49 | whereneg = where(flux_check(lo,la,*) lt 0) ; (12 pts max) Identification de potentiels points à problÚmes, corrigés négativement |
---|
50 | nbannul=n_elements(whereneg)*(total(whereneg) ne -1) |
---|
51 | flux_corr=flux_check(lo,la,*) ; Création d'un vecteur pour recevoir les valeurs corrigées, initialisé à flux_check au cas où on n'ait rien à faire d'autre qu'une seule itération |
---|
52 | A2 = A ; Je repars de la matrice A initiale, ce pour chaque point de grille |
---|
53 | ; ; Potentiellement plusieurs passages pour éliminer toutes les valeurs négatives |
---|
54 | while nbannul ne 0 do begin ; Si l'on a effectivement des émissions corrigées négativement... |
---|
55 | m_bloq(lo,la,whereneg) = 1 ; Update de la matrice m_bloq |
---|
56 | for m=0,11 do begin |
---|
57 | if m eq 11 then begin ; Pour plus de facilité, mois précédents et suivants codés ici |
---|
58 | p=10 |
---|
59 | s=0 |
---|
60 | endif else if m eq 0 then begin |
---|
61 | p=11 |
---|
62 | s=1 |
---|
63 | endif else begin |
---|
64 | p = m-1 |
---|
65 | s = m+1 |
---|
66 | endelse |
---|
67 | if m_bloq(lo,la,m) then begin ; Je traite les mois bloqués en eux-mêmes |
---|
68 | A2(p,m) = 0. |
---|
69 | A2(m,m) = 1. |
---|
70 | A2(s,m) = 0. |
---|
71 | endif ; Fin du cas si l'on est sur un mois bloqué |
---|
72 | if ~m_bloq(lo,la,m) then begin ; Je traite les mois non bloqués, pour ceux adjacents à un mois bloqué |
---|
73 | if m_bloq(lo,la,p) and m_bloq(lo,la,s) then begin ; Mois encadré de deux mois bloqués |
---|
74 | A2(p,m) = 1./4. |
---|
75 | A2(m,m) = 1./2. |
---|
76 | A2(s,m) = 1./4. |
---|
77 | endif else if m_bloq(lo,la,p) and ~m_bloq(lo,la,s) then begin ; Mois précédent bloqué (uniquement) |
---|
78 | A2(p,m) = 2./8. |
---|
79 | A2(m,m) = 5./8. |
---|
80 | A2(s,m) = 1./8. |
---|
81 | endif else if ~m_bloq(lo,la,p) and m_bloq(lo,la,s) then begin ; Mois suivant bloqué (uniquement) |
---|
82 | A2(p,m) = 1./8. |
---|
83 | A2(m,m) = 5./8. |
---|
84 | A2(s,m) = 2./8. |
---|
85 | endif |
---|
86 | endif ; Fin du cas mois non bloqué |
---|
87 | endfor ; Fin de la boucle sur les mois, balayage de la matrice |
---|
88 | flux_corr = invert(A2) ## transpose(fluxinit(lo,la,*)) ; Ré-itération de la multiplication matricielle avec la matrice A modifiée (A2) |
---|
89 | whereneg = where(flux_corr lt 0) ; (12 pts max) Ré-identification de potentiels nouveaux points à problÚmes, corrigés négativement |
---|
90 | nbannul=n_elements(whereneg)*(total(whereneg) ne -1) |
---|
91 | endwhile ; Fin du cas où l'on avait des problÚmes d'émissions corrigées négativement |
---|
92 | ; *** IMPORTANT ! *** Pour signaler les mois bloqués, on prend la convention suivante : |
---|
93 | ; valeur négative ou nulle <=> mois bloqué |
---|
94 | ; valeur positive <=> mois à interpolation classique |
---|
95 | flux(lo,la,*) = flux_corr ; En sortie de la boucle while, normalement flux_corr est complÚtement positif |
---|
96 | endfor ; Fin boucle lat |
---|
97 | endfor ; Fin boucle lon |
---|
98 | nbnegtotal = n_elements(where(m_bloq eq 1)) * (total(where(m_bloq eq 1)) ne -1) |
---|
99 | if nbnegtotal ne 0 then flux(where(m_bloq eq 1)) = -flux(where(m_bloq eq 1)) ; Je force à des valeurs négatives |
---|
100 | ; |
---|
101 | month_in_year=12 |
---|
102 | nbpoint=1682 |
---|
103 | flux2=fltarr(nbpoint,month_in_year) |
---|
104 | flux2(*,*)=0.0 |
---|
105 | ; |
---|
106 | for l=0,month_in_year-1 do begin |
---|
107 | flux2(0,l)=TOTAL(flux(*,0,l))/float(dimlon) |
---|
108 | flux2(nbpoint-1,l)=TOTAL(flux(*,dimlat-1,l))/float(dimlon) |
---|
109 | endfor |
---|
110 | ; |
---|
111 | pt=1 |
---|
112 | for j=1,dimlat-2 do begin |
---|
113 | for i=0,dimlon-1 do begin |
---|
114 | for l=0,month_in_year-1 do begin |
---|
115 | flux2(pt,l)=flux(i,j,l) |
---|
116 | endfor |
---|
117 | pt=pt+1 |
---|
118 | endfor |
---|
119 | endfor |
---|
120 | ; |
---|
121 | ;saving netcdf file |
---|
122 | ; |
---|
123 | fluxstruct={vector:fltarr(nbpoint),time:fltarr(month_in_year), $ |
---|
124 | flx_CO2:fltarr(nbpoint,month_in_year) } |
---|
125 | ; |
---|
126 | fluxstruct.vector=float(indgen(nbpoint)+1) |
---|
127 | ;;fluxstruct.time=float(indgen(month_in_year)+1) |
---|
128 | fluxstruct.time=[15, 45, 75, 105, 135, 165, 195, 225, 255, 285, 315, 345] |
---|
129 | fluxstruct.flx_CO2=flux2 |
---|
130 | ; |
---|
131 | attributes = {units:strarr(3),long_name:strarr(3)} |
---|
132 | attributes.units = ['vector','days since 1960-01-01','flux'] |
---|
133 | attributes.long_name = ['vector', 'time', 'flux'] |
---|
134 | ; |
---|
135 | dimensions = {isdim:intarr(3), links:intarr(2,3)} |
---|
136 | dimensions.isdim = [1,1,0] ; (1=dimension, 0=variable) |
---|
137 | dimensions.links = [ [-1,-1],[-1,-1],[0,1] ] |
---|
138 | ; |
---|
139 | netcdfwrite,'/data/obolmd/CMIP6/CO2_EMISSIONS_48x36/flux_1D_CO2_1867.nc',fluxstruct,clobber=1, attributes=attributes, dimensions=dimensions |
---|
140 | ; |
---|
141 | end |
---|