1 | PRO bining2, $ ; |
---|
2 | ;INPUTS |
---|
3 | density, $ ; density at T points (3D array) |
---|
4 | x1, $ ; field at T points (e.g. Salinity or Temperature) |
---|
5 | sobwlmax, $ ; bowl depth array |
---|
6 | sig_bowl, $ ; switch for bowl overlay |
---|
7 | ;OUTPUTS |
---|
8 | depth_bin, $ ; depth of layers (3D array) |
---|
9 | thick_bin, $ ; thickness of layers (3D array) |
---|
10 | x1_bin, $ ; X averaged for each sigma-layer (3D array) |
---|
11 | bowl_bin, $ ; bowl depth binned on density |
---|
12 | ; KEYWORDS |
---|
13 | SIGMA = sigma, $ ; bining values (optionnal) |
---|
14 | DEPTH_T = depth_t, $ ; depth of T level (NECESSARY) |
---|
15 | DEPTH_W = depth_w, $ ; depth of W level (NECESSARY) |
---|
16 | E3T = e3t, $ |
---|
17 | E3W = e3w, $ |
---|
18 | TMASK = tmask ; tmask (3D array) (NECESSARY) |
---|
19 | |
---|
20 | ;; |
---|
21 | ;; PRINCIPE DU BINING |
---|
22 | ;; |
---|
23 | ;; - etant donnee une liste de densite (sigma) on cherche les profondeurs |
---|
24 | ;; correspondants a ces valeurs. Pour cela on utilise une interpolation lineaire. |
---|
25 | ;; |
---|
26 | ;; - puis, connaissant les profondeurs des isopycnes, on interpole les valeurs du contenu |
---|
27 | ;; vertical de x a ces profondeurs. Il est preferable d' interpoler le contenu de x plutot |
---|
28 | ;; que x afin d'assurer la conservation de x. |
---|
29 | ;; |
---|
30 | ;; - on deduit la valeur moyenne de x dans chaque couche (i.e. entre deux isopycnes). Il ne |
---|
31 | ;; s'agit donc pas de la projection de x selon les isopycnes mais de la valeur moyenne entre |
---|
32 | ;; deux isopycnes. Les proprietes sont meilleures (conservation de la quantite totale de x). |
---|
33 | ;; |
---|
34 | ;; REMARQUES |
---|
35 | ;; |
---|
36 | ;; - il n'y a aucun appel a un fichier de common (d'ou la liste des mots clef un peu longue) |
---|
37 | ;; |
---|
38 | ;; - le calcul n'est rigoureux que pour les points T. Le bining de la vitesse doit etre adapte... |
---|
39 | ;; |
---|
40 | ;; - si l'utilisateur donne N interfaces, cela defini N+1 couches, d'ou les decalages dans les indices |
---|
41 | ;; et l'apparition d'un +1 dans la declaration des tableaux de sortie |
---|
42 | ;; la couche 0 contient tous les points dont la densite est inferieure a la densite minimale du bining |
---|
43 | ;; la couche N+1 " " " superieure " maximale |
---|
44 | ;; |
---|
45 | ;; |
---|
46 | ;; Creation : 19/10/99 G. Roullet |
---|
47 | ;; |
---|
48 | size3d = size(density) |
---|
49 | jpi = size3d(1) |
---|
50 | jpj = size3d(2) |
---|
51 | jpk = size3d(3) |
---|
52 | ;print, 'size3d=', size3d |
---|
53 | ;print, 'size(tmask)=', size(tmask) |
---|
54 | |
---|
55 | indm = where(tmask eq 0) |
---|
56 | |
---|
57 | x1(indm) = !VALUES.F_NAN |
---|
58 | density(indm) = !VALUES.F_NAN |
---|
59 | |
---|
60 | N_z = jpk |
---|
61 | |
---|
62 | ; adjust profiles to avoid static instabilities (there are almost none IF using pac_40) |
---|
63 | IF jpi EQ 180 THEN BEGIN |
---|
64 | density = npc(density) |
---|
65 | ENDIF |
---|
66 | |
---|
67 | ; N_s, defini ci dessous est le nombre d'interface |
---|
68 | |
---|
69 | IF keyword_set(sigma) THEN BEGIN |
---|
70 | s_s = sigma |
---|
71 | N_s = n_elements(s_s) |
---|
72 | ENDIF ELSE BEGIN |
---|
73 | ; |
---|
74 | ; Definition d'un bining arbitraire (si density est un sigma potentiel, cela repere |
---|
75 | ; les eaux intermediaires) |
---|
76 | ; |
---|
77 | N_s = 20 |
---|
78 | s_s = 27.6+0.2*findgen(N_s)/(N_s-1) |
---|
79 | N_s = 3 |
---|
80 | s_s = [26.8, 27.6, 27.8] |
---|
81 | s_s = [22.8, 27.6, 27.8] |
---|
82 | ENDELSE |
---|
83 | |
---|
84 | ; print, 'Density bining :', s_s |
---|
85 | ;; |
---|
86 | ;; Definition des variables |
---|
87 | ;; |
---|
88 | ; Profils selon les niveaux du modele (suffixe _z) |
---|
89 | ; |
---|
90 | c1_z = fltarr(N_z) ; profil du contenu vertical de x |
---|
91 | s_z = fltarr(N_z) ; profil de la densite |
---|
92 | z_zt = depth_t ; profondeur au point T (k=0 -> 5m) |
---|
93 | z_zw = depth_w ; W (k=0 -> 10m) |
---|
94 | |
---|
95 | ; Profils selon les couches de densite (suffixe _s) |
---|
96 | ; |
---|
97 | z_s = fltarr(N_s+1) |
---|
98 | c1_s = fltarr(N_s+1) |
---|
99 | x1_s = fltarr(N_s+1) |
---|
100 | bowl_s = 0. |
---|
101 | |
---|
102 | ; Tableaux de sorties |
---|
103 | ; |
---|
104 | depth_bin = fltarr(jpi, jpj, N_s+1) |
---|
105 | thick_bin = fltarr(jpi, jpj, N_s+1) |
---|
106 | x1_bin = fltarr(jpi, jpj, N_s+1) |
---|
107 | bowl_bin = fltarr(jpi, jpj) |
---|
108 | |
---|
109 | x1_bin(*, *, *) = !VALUES.F_NAN |
---|
110 | depth_bin(*, *, *) = !VALUES.F_NAN |
---|
111 | thick_bin(*, *, *) = !VALUES.F_NAN |
---|
112 | bowl_bin(*, *) = !VALUES.F_NAN |
---|
113 | |
---|
114 | |
---|
115 | ; Calcul du contenu de x sur la verticale (permet d assurer la conservation integrale) |
---|
116 | ; |
---|
117 | x1_content = fltarr(jpi, jpj, jpk) |
---|
118 | |
---|
119 | |
---|
120 | ; x1_content(*, *, 0) = e3t(0) * x1(*, *, 0) * tmask (*, *, 0) |
---|
121 | ; FOR k = 1, jpk-1 DO x1_content(*, *, k) = x1_content(*, *, k-1) $ |
---|
122 | ; + e3t(k) * x1(*, *, k)* tmask (*, *, k) |
---|
123 | x1_content = x1 |
---|
124 | |
---|
125 | ; |
---|
126 | ; Boucle sur le domaine 2D |
---|
127 | ; |
---|
128 | |
---|
129 | FOR i = 0, (jpi-1) DO BEGIN |
---|
130 | ; FOR i = 20, 20 DO BEGIN |
---|
131 | ; FOR i = 162, 162 DO BEGIN |
---|
132 | FOR j = 0, (jpj-1) DO BEGIN |
---|
133 | ; FOR j = 20, 20 DO BEGIN |
---|
134 | ; FOR j = 116, 116 DO BEGIN |
---|
135 | ; |
---|
136 | ; Indices des points T dans l ocean |
---|
137 | ;print, ' ' |
---|
138 | ;print, ' ===============================' |
---|
139 | ;print, ' POINT I,J =', i, j |
---|
140 | ; |
---|
141 | ;stop |
---|
142 | i_ocean = where(tmask(i, j, *) EQ 1) |
---|
143 | ;print, 'tmask(i, j, *)=', tmask(i, j, *) |
---|
144 | ;print, 'i_ocean=', i_ocean |
---|
145 | z_s = z_s*0. |
---|
146 | c1_s(*) = !VALUES.F_NAN |
---|
147 | bowl_s = !VALUES.F_NAN |
---|
148 | ; x1_s(*) = !VALUES.F_NAN |
---|
149 | |
---|
150 | |
---|
151 | IF i_ocean[0] NE -1 THEN BEGIN ; on n entre que si il y a des points ocean |
---|
152 | ; print, 'ocean point' |
---|
153 | i_bottom = i_ocean(n_elements(i_ocean)-1) |
---|
154 | ;print, 'i_bottom=', i_bottom |
---|
155 | z_s(N_s) = z_zw(i_bottom) |
---|
156 | c1_s(N_s) = x1_content(i, j, jpk-1) |
---|
157 | |
---|
158 | s_z(*) = density(i, j, *) |
---|
159 | c1_z(*) = x1_content(i, j, *) |
---|
160 | ; print, 'density profile s_z', s_z |
---|
161 | ; print, 'field profile c1_z', c1_z |
---|
162 | |
---|
163 | ; print, 'bins s_s', s_s |
---|
164 | |
---|
165 | ; extraction d'un sous profil strictement croissant |
---|
166 | ; print, 's_z(i_ocean)=', s_z(i_ocean) |
---|
167 | mini = min(s_z(i_ocean)) |
---|
168 | maxi = max(s_z(i_ocean)) |
---|
169 | i_min = where(s_z(i_ocean) EQ mini) |
---|
170 | i_max = where(s_z(i_ocean) EQ maxi) |
---|
171 | ; print, 'mini, maxi', mini, maxi |
---|
172 | ; print, 'i_min, i_max', i_min, i_max |
---|
173 | ; on prend le plus grand des indices min |
---|
174 | ; et le plus petit des indices max |
---|
175 | i_min = i_min(n_elements(i_min)-1) |
---|
176 | i_max = i_max[0] |
---|
177 | ; print, 'i_min, i_max', i_min, i_max |
---|
178 | ;IF i_min GT i_max THEN print, 'WARNING: i_min, i_max=', i_min, i_max, ' at i,j=,', i, j |
---|
179 | ; on prend le plus grand des indices min |
---|
180 | |
---|
181 | ; IF i_max GE jpk-1 THEN print, i, j, i_max |
---|
182 | ; IF i_min LE 1 THEN print, i, j, i_min |
---|
183 | |
---|
184 | ; Si la valeur du niveau (s_s) est plus faible que la densite de surface, |
---|
185 | ; l isopycne est mise en surface (z_s=0) |
---|
186 | ; |
---|
187 | ; print, 's_z(i_min)', s_z(i_min) |
---|
188 | ;print, 's_s=', s_s |
---|
189 | ind = where(s_s LT s_z(i_min)) |
---|
190 | IF ind[0] NE -1 THEN BEGIN |
---|
191 | ; IF i_min GT i_max THEN print, 'min reached at sigma indices', ind |
---|
192 | z_s(ind) = 0 |
---|
193 | c1_s(ind) = !VALUES.F_NAN |
---|
194 | ENDIF |
---|
195 | |
---|
196 | ; Si la valeur du niveau (s_s) est plus elevee que la densite du fond, |
---|
197 | ; l isopycne est mise au fond (z_s=z_zw(i_bottom)) |
---|
198 | ; |
---|
199 | ; print, 's_z(i_max)', s_z(i_max) |
---|
200 | ind = where(s_s GT s_z(i_max)) |
---|
201 | |
---|
202 | IF ind[0] NE -1 THEN BEGIN |
---|
203 | ; IF i_min GT i_max THEN print, 'max reached at sigma indices', ind |
---|
204 | z_s(ind) = z_s(N_s) |
---|
205 | c1_s(ind) = c1_s(N_s) |
---|
206 | c1_s(ind) = !VALUES.F_NAN |
---|
207 | ENDIF |
---|
208 | ; cas general |
---|
209 | ind = where( (s_s GE s_z(i_min)) AND (s_s LE s_z(i_max)) ) |
---|
210 | ; IF i_min GT i_max THEN print, 's_s indices inside min/max', ind |
---|
211 | IF ind[0] NE -1 THEN BEGIN |
---|
212 | i_profil = i_ocean(i_min:i_max) ; problem line |
---|
213 | |
---|
214 | z_s(ind) = interpol(z_zt(i_profil), s_z(i_profil), s_s(ind)) ; original |
---|
215 | ; z_s(ind) = interpol(z_zw(i_profil), s_z(i_profil), s_s(ind)) ; changed to z_zw 1/7/04 |
---|
216 | |
---|
217 | ; |
---|
218 | ; j'utilise la fonction spline pour interpoler le contenu |
---|
219 | ; l'interpolation lineaire marche aussi. Elle donne des resultats differents. Elle |
---|
220 | ; me semble moins propre. Je la donne en remarque. |
---|
221 | ; |
---|
222 | IF n_elements(i_profil) GT 2 THEN BEGIN |
---|
223 | ; c1_s(ind) = spline(z_zw(i_profil), c1_z(i_profil), z_s(ind), 1) ; cubic spline |
---|
224 | c1_s(ind) = interpol(c1_z(i_profil), z_zt(i_profil), z_s(ind)) ; linear interpolation |
---|
225 | ; SPLINE_P, z_zw(i_profil), c1_z(i_profil), z_s(ind), c1_s(ind) ; generalization of spline(*) |
---|
226 | ;print, z_zw(i_profil), c1_z(i_profil), z_s(ind), c1_s(ind) |
---|
227 | ENDIF ELSE BEGIN |
---|
228 | c1_s(ind) = interpol(c1_z(i_profil), z_zt(i_profil), z_s(ind)) |
---|
229 | ENDELSE |
---|
230 | ; |
---|
231 | ; bowl depth binning |
---|
232 | |
---|
233 | IF sig_bowl EQ 1 THEN BEGIN |
---|
234 | bowl_s = interpol(s_z(i_profil), z_zt(i_profil), sobwlmax(i, j)) |
---|
235 | ENDIF ELSE bowl_s = 0 |
---|
236 | |
---|
237 | ;stop |
---|
238 | ; |
---|
239 | ; |
---|
240 | ; x1_s(ind) = (c1_s(ind+1)-c1_s(ind))/(z_s(ind+1)-z_s(ind)) |
---|
241 | ; x1_s = c1_s |
---|
242 | |
---|
243 | ENDIF |
---|
244 | ;ELSE print, 'ind[0] = -1', ind, i_bottom |
---|
245 | ENDIF |
---|
246 | ;ELSE print, ' land ', tmask(i, j, 1) |
---|
247 | depth_bin (i, j, *) = z_s |
---|
248 | thick_bin (i, j, 0) = z_s(0) |
---|
249 | thick_bin (i, j, 1:N_s) = z_s(1:N_s)-z_s(0:N_s-1) |
---|
250 | IF total(thick_bin (i, j, 1:N_s) LT 0) GT 0 THEN BEGIN |
---|
251 | ;print, 'WARNING: negative thick_bin (i, j, 1:N_s) at i,j= ', i, j |
---|
252 | ;print, 'depth_bin (i, j, 1:N_s)= ', depth_bin (i, j, 1:N_s) |
---|
253 | ;print, 'thick_bin (i, j, 1:N_s)= ', thick_bin (i, j, 1:N_s) |
---|
254 | ;print, 's_z= ', s_z |
---|
255 | ;stop |
---|
256 | ENDIF |
---|
257 | x1_bin (i, j, *) = c1_s |
---|
258 | bowl_bin (i, j) = bowl_s |
---|
259 | ; print, 'c1_s', c1_s |
---|
260 | ; |
---|
261 | ENDFOR |
---|
262 | ENDFOR |
---|
263 | ; mask depth_bin with undefined values of x_bin |
---|
264 | depth_bin(where(finite(x1_bin, /nan) EQ 1)) = !VALUES.F_NAN |
---|
265 | thick_bin(where(finite(x1_bin, /nan) EQ 1)) = !VALUES.F_NAN |
---|
266 | |
---|
267 | ; OPTIONAL: compute spiciness by removing isopycnal mean for each sigma (T or S) |
---|
268 | IF 0 THEN BEGIN |
---|
269 | print, ' Computing spiciness...' |
---|
270 | FOR i = 0, (size(x1_bin))[3]-1 DO BEGIN |
---|
271 | x1_bin(*,*,i) = x1_bin(*,*,i) - mean(x1_bin(*,*,i),/NaN) |
---|
272 | ENDFOR |
---|
273 | ENDIF |
---|
274 | |
---|
275 | ; OPTIONAL: remove domain mean (S) |
---|
276 | IF 0 THEN BEGIN |
---|
277 | print, ' Removing domain mean...' |
---|
278 | x1_bin(*,*,*) = x1_bin(*,*,*) - mean(x1_bin(*,*,*),/NaN) |
---|
279 | ENDIF |
---|
280 | |
---|
281 | ; OPTIONAL: remove 34.6psu (S) |
---|
282 | IF 1 THEN BEGIN |
---|
283 | print, ' Removing 34.6psu...' |
---|
284 | x1_bin(*,*,*) = x1_bin(*,*,*) - 34.6 |
---|
285 | ENDIF |
---|
286 | |
---|
287 | END |
---|