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 | ;OUTPUTS |
---|
6 | depth_bin, $ ; depth of layers (3D array) |
---|
7 | thick_bin, $ ; thickness of layers (3D array) |
---|
8 | x1_bin, $ ; X averaged for each sigma-layer (3D array) |
---|
9 | ; KEYWORDS |
---|
10 | SIGMA = sigma, $ ; bining values (optionnal) |
---|
11 | DEPTH_T = depth_t, $ ; depth of T level (NECESSARY) |
---|
12 | DEPTH_W = depth_w, $ ; depth of W level (NECESSARY) |
---|
13 | E3T = e3t, $ |
---|
14 | E3W = e3w, $ |
---|
15 | TMASK = tmask ; tmask (3D array) (NECESSARY) |
---|
16 | |
---|
17 | ;; |
---|
18 | ;; PRINCIPE DU BINING |
---|
19 | ;; |
---|
20 | ;; - etant donnee une liste de densite (sigma) on cherche les profondeurs |
---|
21 | ;; correspondants a ces valeurs. Pour cela on utilise une interpolation lineaire. |
---|
22 | ;; |
---|
23 | ;; - puis, connaissant les profondeurs des isopycnes, on interpole les valeurs du contenu |
---|
24 | ;; vertical de x a ces profondeurs. Il est preferable d' interpoler le contenu de x plutot |
---|
25 | ;; que x afin d'assurer la conservation de x. |
---|
26 | ;; |
---|
27 | ;; - on deduit la valeur moyenne de x dans chaque couche (i.e. entre deux isopycnes). Il ne |
---|
28 | ;; s'agit donc pas de la projection de x selon les isopycnes mais de la valeur moyenne entre |
---|
29 | ;; deux isopycnes. Les proprietes sont meilleures (conservation de la quantite totale de x). |
---|
30 | ;; |
---|
31 | ;; REMARQUES |
---|
32 | ;; |
---|
33 | ;; - il n'y a aucun appel a un fichier de common (d'ou la liste des mots clef un peu longue) |
---|
34 | ;; |
---|
35 | ;; - le calcul n'est rigoureux que pour les points T. Le bining de la vitesse doit etre adapte... |
---|
36 | ;; |
---|
37 | ;; - si l'utilisateur donne N interfaces, cela defini N+1 couches, d'ou les decalages dans les indices |
---|
38 | ;; et l'apparition d'un +1 dans la declaration des tableaux de sortie |
---|
39 | ;; la couche 0 contient tous les points dont la densite est inferieure a la densite minimale du bining |
---|
40 | ;; la couche N+1 " " " superieure " maximale |
---|
41 | ;; |
---|
42 | ;; |
---|
43 | ;; Creation : 19/10/99 G. Roullet |
---|
44 | ;; |
---|
45 | size3d = size(density) |
---|
46 | jpi = size3d(1) |
---|
47 | jpj = size3d(2) |
---|
48 | jpk = size3d(3) |
---|
49 | |
---|
50 | indm = where(tmask eq 0) |
---|
51 | |
---|
52 | x1(indm) = 0. |
---|
53 | density(indm) = !VALUES.F_NAN |
---|
54 | |
---|
55 | N_z = jpk |
---|
56 | |
---|
57 | ; |
---|
58 | ; N_s, defini ci dessous est le nombre d'interface |
---|
59 | |
---|
60 | IF keyword_set(sigma) THEN BEGIN |
---|
61 | s_s = sigma |
---|
62 | N_s = n_elements(s_s) |
---|
63 | ENDIF ELSE BEGIN |
---|
64 | ; |
---|
65 | ; Definition d'un bining arbitraire (si density est un sigma potentiel, cela repere |
---|
66 | ; les eaux intermediaires) |
---|
67 | ; |
---|
68 | N_s = 20 |
---|
69 | s_s = 27.6+0.2*findgen(N_s)/(N_s-1) |
---|
70 | N_s = 3 |
---|
71 | s_s = [26.8, 27.6, 27.8] |
---|
72 | s_s = [22.8, 27.6, 27.8] |
---|
73 | ENDELSE |
---|
74 | |
---|
75 | ; print, 'Density bining :', s_s |
---|
76 | ;; |
---|
77 | ;; Definition des variables |
---|
78 | ;; |
---|
79 | ; Profils selon les niveaux du modele (suffixe _z) |
---|
80 | ; |
---|
81 | c1_z = fltarr(N_z) ; profil du contenu vertical de x |
---|
82 | s_z = fltarr(N_z) ; profil de la densite |
---|
83 | z_zt = depth_t ; profondeur au point T (k=0 -> 5m) |
---|
84 | z_zw = depth_w ; W (k=0 -> 10m) |
---|
85 | |
---|
86 | ; Profils selon les couches de densite (suffixe _s) |
---|
87 | ; |
---|
88 | z_s = fltarr(N_s+1) |
---|
89 | c1_s = fltarr(N_s+1) |
---|
90 | x1_s = fltarr(N_s+1) |
---|
91 | |
---|
92 | ; Tableaux de sorties |
---|
93 | ; |
---|
94 | depth_bin = fltarr(jpi, jpj, N_s+1) |
---|
95 | thick_bin = fltarr(jpi, jpj, N_s+1) |
---|
96 | x1_bin = fltarr(jpi, jpj, N_s+1) |
---|
97 | |
---|
98 | x1_bin(*, *, *) = !VALUES.F_NAN |
---|
99 | depth_bin(*, *, *) = !VALUES.F_NAN |
---|
100 | thick_bin(*, *, *) = !VALUES.F_NAN |
---|
101 | |
---|
102 | ; |
---|
103 | ; Calcul du contenu de x sur la verticale (permet d assurer la conservation integrale) |
---|
104 | ; |
---|
105 | x1_content = fltarr(jpi, jpj, jpk) |
---|
106 | x1_content(*, *, 0) = e3t(0) * x1(*, *, 0) * tmask (*, *, 0) |
---|
107 | FOR k = 1, jpk-1 DO x1_content(*, *, k) = x1_content(*, *, k-1) $ |
---|
108 | + e3t(k) * x1(*, *, k)* tmask (*, *, k) |
---|
109 | |
---|
110 | |
---|
111 | ; |
---|
112 | ; Boucle sur le domaine 2D |
---|
113 | ; |
---|
114 | FOR i = 0, (jpi-1) DO BEGIN |
---|
115 | FOR j = 0, (jpj-1) DO BEGIN |
---|
116 | ; |
---|
117 | ; Indices des points T dans l ocean |
---|
118 | ; |
---|
119 | i_ocean = where(tmask(i, j, *) EQ 1) |
---|
120 | |
---|
121 | z_s = z_s*0. |
---|
122 | c1_s = c1_s*0. |
---|
123 | x1_s = x1_s*0. |
---|
124 | |
---|
125 | |
---|
126 | IF i_ocean[0] NE -1 THEN BEGIN ; on n entre que si il y a des points ocean |
---|
127 | |
---|
128 | i_bottom = i_ocean(n_elements(i_ocean)-1) |
---|
129 | |
---|
130 | z_s(N_s) = z_zw(i_bottom) |
---|
131 | c1_s(N_s) = x1_content(i, j, jpk-1) |
---|
132 | |
---|
133 | s_z(*) = density(i, j, *) |
---|
134 | c1_z(*) = x1_content(i, j, *) |
---|
135 | |
---|
136 | ; extraction d'un sous profil strictement croissant |
---|
137 | mini = min(s_z(i_ocean)) |
---|
138 | maxi = max(s_z(i_ocean)) |
---|
139 | i_min = where(s_z(i_ocean) EQ mini) |
---|
140 | i_max = where(s_z(i_ocean) EQ maxi) |
---|
141 | ; on prend le plus grand des indices min |
---|
142 | ; et le plus petit des indices max |
---|
143 | i_min = i_min(n_elements(i_min)-1) |
---|
144 | i_max = i_max[0] |
---|
145 | |
---|
146 | ; IF i_max GE jpk-1 THEN print, i, j, i_max |
---|
147 | ; IF i_min LE 1 THEN print, i, j, i_min |
---|
148 | |
---|
149 | ; Si la valeur du niveau (s_s) est plus faible que la densite de surface, |
---|
150 | ; l isopycne est mise en surface (z_s=0) |
---|
151 | ; |
---|
152 | ind = where(s_s LT s_z(i_min)) |
---|
153 | IF ind[0] NE -1 THEN BEGIN |
---|
154 | z_s(ind) = 0 |
---|
155 | c1_s(ind) = 0 |
---|
156 | ENDIF |
---|
157 | |
---|
158 | ; Si la valeur du niveau (s_s) est plus elevee que la densite du fond, |
---|
159 | ; l isopycne est mise au fond (z_s=z_zw(i_bottom)) |
---|
160 | ; |
---|
161 | ind = where(s_s GT s_z(i_max)) |
---|
162 | IF ind[0] NE -1 THEN BEGIN |
---|
163 | z_s(ind) = z_s(N_s) |
---|
164 | c1_s(ind) = c1_s(N_s) |
---|
165 | ENDIF |
---|
166 | ; cas general |
---|
167 | ind = where( (s_s GE s_z(i_min)) AND (s_s LE s_z(i_max)) ) |
---|
168 | IF ind[0] NE -1 THEN BEGIN |
---|
169 | |
---|
170 | i_profil = i_ocean(i_min:i_max) |
---|
171 | |
---|
172 | z_s(ind) = interpol(z_zt(i_profil), s_z(i_profil), s_s(ind)) |
---|
173 | |
---|
174 | ; |
---|
175 | ; j'utilise la fonction spline pour interpoler le contenu |
---|
176 | ; l'interpolation lineaire marche aussi. Elle donne des resultats differents. Elle |
---|
177 | ; me semble moins propre. Je la donne en remarque. |
---|
178 | ; |
---|
179 | IF n_elements(i_profil) GT 2 THEN BEGIN |
---|
180 | c1_s(ind) = spline(z_zw(i_profil), c1_z(i_profil), z_s(ind), 1) |
---|
181 | ; c1_s(ind) = interpol(c1_z(i_profil), z_zw(i_profil), z_s(ind)) |
---|
182 | ENDIF ELSE BEGIN |
---|
183 | c1_s(ind) = interpol(c1_z(i_profil), z_zw(i_profil), z_s(ind)) |
---|
184 | ENDELSE |
---|
185 | ; |
---|
186 | ; |
---|
187 | x1_s(ind) = (c1_s(ind+1)-c1_s(ind))/(z_s(ind+1)-z_s(ind)) |
---|
188 | |
---|
189 | ENDIF |
---|
190 | ENDIF |
---|
191 | ; |
---|
192 | depth_bin (i, j, *) = z_s |
---|
193 | thick_bin (i, j, 0) = z_s(0) |
---|
194 | thick_bin (i, j, 1:N_s) = z_s(1:N_s)-z_s(0:N_s-1) |
---|
195 | x1_bin (i, j, *) = x1_s |
---|
196 | ; |
---|
197 | ENDFOR |
---|
198 | ENDFOR |
---|
199 | ; mask depth_bin with undefined values of x_bin |
---|
200 | depth_bin(where(finite(x1_bin, /nan) EQ 1)) = !VALUES.F_NAN |
---|
201 | thick_bin(where(finite(x1_bin, /nan) EQ 1)) = !VALUES.F_NAN |
---|
202 | |
---|
203 | |
---|
204 | END |
---|