[2] | 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 |
---|