;+ ; PRINCIPE DU BINING ; ; - etant donnee une liste de densite (sigma) on cherche les profondeurs ; correspondants a ces valeurs. Pour cela on utilise une interpolation lineaire. ; ; - puis, connaissant les profondeurs des isopycnes, on interpole les valeurs du contenu ; vertical de x a ces profondeurs. Il est preferable d' interpoler le contenu de x plutot ; que x afin d'assurer la conservation de x. ; ; - on deduit la valeur moyenne de x dans chaque couche (i.e. entre deux isopycnes). Il ne ; s'agit donc pas de la projection de x selon les isopycnes mais de la valeur moyenne entre ; deux isopycnes. Les proprietes sont meilleures (conservation de la quantite totale de x). ; ; REMARQUES ; ; - il n'y a aucun appel a un fichier de common (d'ou la liste des mots clef un peu longue) ; ; - le calcul n'est rigoureux que pour les points T. Le bining de la vitesse doit etre adapte... ; ; - si l'utilisateur donne N interfaces, cela defini N+2 couches, d'ou les decalages dans les indices ; et l'apparition d'un +2 dans la declaration des tableaux de sortie ; la couche 0 contient tous les points dont la densite est inferieure a la densite minimale du bining ; la couche N+1 " " " superieure " maximale ; ; @param DENSITY {in} ; density at T points (3D array) ; ; @param X {in} ; field at T points (e.g. Salinity or Temperature) ; ; @param DEPTH_BIN {out} ; depth of sopycnical surfaces (3D array) ; ; @param X_BIN {out} ; X averaged for each sigma-layer (3D array) ; ; @keyword SIGMA {in}{optional} ; bining values ; ; @keyword DEPTH_T {in}{required} ; depth of T level ; ; @keyword DEPTH_W {in}{required} ; depth of W level ; ; @keyword E3T {in}{required} ; ; @keyword TMASK {in}{required} ; tmask (3D array) ; ; @history ; 19/10/99 G. Roullet ; ; @version ; $Id$ ; ;- PRO bining $ , density $ , x $ , depth_bin $ , x_bin $ , SIGMA=sigma $ , DEPTH_T=depth_t $ , DEPTH_W=depth_w $ , E3T=e3t $ , TMASK=tmask ; compile_opt idl2, strictarrsubs ; size3d = size(density) jpi = size3d[1] jpj = size3d[2] jpk = size3d[3] indm = where(tmask eq 0) x[indm] = 0. density[indm] = !VALUES.F_NAN N_z = jpk ; ; N_s, defini ci dessous est le nombre d'interface IF keyword_set(sigma) THEN BEGIN s_s = sigma N_s = n_elements(s_s) ENDIF ELSE BEGIN ; ; Definition d'un bining arbitraire (si density est un sigma potentiel, cela repere ; les eaux intermediaires) ; N_s = 20 s_s = 27.6+0.2*findgen(N_s)/(N_s-1) ENDELSE ;; ;; Definition des variables ;; ; Profils selon les niveaux du modele (suffixe _z) ; c_z = fltarr(N_z) ; profil du contenu vertical de x s_z = fltarr(N_z) ; profil de la densite z_zt = depth_t ; profondeur au point T (k=0 -> 5m) z_zw = depth_w ; W (k=0 -> 10m) ; Profils selon les couches de densite (suffixe _s) ; z_s = fltarr(N_s+2) c_s = fltarr(N_s+2) x_s = fltarr(N_s+2) ; Tableaux de sorties ; depth_bin = fltarr(jpi, jpj, N_s+2) x_bin = fltarr(jpi, jpj, N_s+2) cont_s = fltarr(jpi, jpj, N_s+2) x_bin[*, *, *] = !VALUES.F_NAN depth_bin[*, *, *] = !VALUES.F_NAN ; ; Calcul du contenu de x sur la verticale (permet d assurer la conservation integrale) ; x_content = fltarr(jpi, jpj, jpk) x_content[*, *, 0] = 0. ; FOR k = 0, jpk-1 DO BEGIN ; print, k, z_zw[k], e3t[k], x[4, 30, k] ; ENDFOR FOR k = 1, jpk-1 DO BEGIN x_content[*, *, k] = x_content[*, *, k-1] $ ; + e3t[k-1] * x[*, *, k-1] + [ z_zw[k]-z_zw[k-1] ] * x[*, *, k] ; print, k, z_zw[k], e3t[k-1], [ z_zw[k]-z_zw[k-1] ], [ z_zw[k]-z_zw[k-1] ]-e3t[k-1] ENDFOR ; ; ; Boucle sur le domaine 2D ; FOR i = 0, (jpi-1) DO BEGIN FOR j = 0, (jpj-1) DO BEGIN ; ; Indices des points T dans l ocean ; i_ocean = where(tmask[i, j, *] EQ 1) z_s[*] = !VALUES.F_NAN c_s = c_s*0. x_s[*] = !VALUES.F_NAN IF i_ocean[0] NE -1 THEN BEGIN ; on n entre que si il y a des points ocean i_bottom = i_ocean[n_elements(i_ocean)-1] z_s[0] = 0 z_s[N_s+1] = z_zw[i_bottom] c_s[0] = 0 c_s[N_s+1] = x_content[i, j, jpk-1] s_z[*] = density[i, j, *] c_z[*] = x_content[i, j, *] ; exctraction d'un sous profil strictement croissant mini = min(s_z[i_ocean]) maxi = max(s_z[i_ocean]) i_min = where(s_z[i_ocean] EQ mini) i_max = where(s_z[i_ocean] EQ maxi) ; on prend le plus grand des indices min ; et le plus petit des indices max i_min = i_min[n_elements(i_min)-1] i_max = i_max[0] ; Si la valeur du niveau (s_s) est plus faible que la densite de surface, ; l isopycne est mise en surface (z_s=0) ; ind = where(s_s LT s_z[i_min]) IF ind[0] NE -1 THEN BEGIN z_s[ind+1] = 0 c_s[ind+1] = 0 ENDIF ; Si la valeur du niveau (s_s) est plus elevee que la densite du fond, ; l isopycne est mise au fond (z_s=z_zw[i_bottom]) ; ind = where(s_s GT s_z[i_max]) IF ind[0] NE -1 THEN BEGIN z_s[ind+1] = z_s[N_s+1] c_s[ind+1] = c_s[N_s+1] ENDIF ind = where( (s_s GE s_z[i_min]) AND (s_s LE s_z[i_max]) ) IF ind[0] NE -1 THEN BEGIN i_profil = i_ocean[i_min:i_max] ; print, i, j, n_elements(i_profil) z_s[ind+1] = interpol(z_zt[i_profil], s_z[i_profil], s_s[ind]) ; ; j'utilise la fonction spline pour interpoler le contenu ; IF n_elements(i_profil) GT 2 THEN BEGIN c_s[ind+1] = spline(z_zw[i_profil], c_z[i_profil], z_s[ind+1], 1) ENDIF ELSE BEGIN c_s[ind+1] = interpol(c_z[i_profil], z_zw[i_profil], z_s[ind+1]) ENDELSE ; ; l'interpolation lineaire marche aussi. Elle donne des resultats differents. Elle ; me semble moins propre. Je la donne en remarque. ; ; c_s[ind+1] = interpol(c_z[i_profil], z_zw[i_profil], z_s[ind+1]) x_s[ind+1] = (c_s[ind+2]-c_s[ind+1])/(z_s[ind+2]-z_s[ind+1]) ENDIF ; x_s[0] = c_s[1]/z_s[1] ENDIF depth_bin[i, j, *] = z_s x_bin [i, j, *] = x_s ENDFOR ENDFOR ; mask depth_bin with undefined values of x_bin depth_bin[where(finite(x_bin, /nan) EQ 1)] = !VALUES.F_NAN ; print, min(x_bin(where (abs(x_bin) LT 1.e19))), max(x_bin(where (abs(x_bin) LT 1.e19))) ; print, min(depth_bin(where (abs(depth_bin) LT 1.e19))), max(depth_bin(where (abs(depth_bin) LT 1.e19))) END