PRO bining2, $ ; ;INPUTS density, $ ; density at T points (3D array) x1, $ ; field at T points (e.g. Salinity or Temperature) ;OUTPUTS depth_bin, $ ; depth of layers (3D array) thick_bin, $ ; thickness of layers (3D array) x1_bin, $ ; X averaged for each sigma-layer (3D array) ; KEYWORDS SIGMA = sigma, $ ; bining values (optionnal) DEPTH_T = depth_t, $ ; depth of T level (NECESSARY) DEPTH_W = depth_w, $ ; depth of W level (NECESSARY) E3T = e3t, $ E3W = e3w, $ TMASK = tmask ; tmask (3D array) (NECESSARY) ;; ;; 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+1 couches, d'ou les decalages dans les indices ;; et l'apparition d'un +1 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 ;; ;; ;; Creation : 19/10/99 G. Roullet ;; size3d = size(density) jpi = size3d(1) jpj = size3d(2) jpk = size3d(3) indm = where(tmask eq 0) x1(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) N_s = 3 s_s = [26.8, 27.6, 27.8] s_s = [22.8, 27.6, 27.8] ENDELSE ; print, 'Density bining :', s_s ;; ;; Definition des variables ;; ; Profils selon les niveaux du modele (suffixe _z) ; c1_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+1) c1_s = fltarr(N_s+1) x1_s = fltarr(N_s+1) ; Tableaux de sorties ; depth_bin = fltarr(jpi, jpj, N_s+1) thick_bin = fltarr(jpi, jpj, N_s+1) x1_bin = fltarr(jpi, jpj, N_s+1) x1_bin(*, *, *) = !VALUES.F_NAN depth_bin(*, *, *) = !VALUES.F_NAN thick_bin(*, *, *) = !VALUES.F_NAN ; ; Calcul du contenu de x sur la verticale (permet d assurer la conservation integrale) ; x1_content = fltarr(jpi, jpj, jpk) x1_content(*, *, 0) = e3t(0) * x1(*, *, 0) * tmask (*, *, 0) FOR k = 1, jpk-1 DO x1_content(*, *, k) = x1_content(*, *, k-1) $ + e3t(k) * x1(*, *, k)* tmask (*, *, k) ; ; 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 = z_s*0. c1_s = c1_s*0. x1_s = x1_s*0. 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(N_s) = z_zw(i_bottom) c1_s(N_s) = x1_content(i, j, jpk-1) s_z(*) = density(i, j, *) c1_z(*) = x1_content(i, j, *) ; extraction 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] ; IF i_max GE jpk-1 THEN print, i, j, i_max ; IF i_min LE 1 THEN print, i, j, i_min ; 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) = 0 c1_s(ind) = 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) = z_s(N_s) c1_s(ind) = c1_s(N_s) ENDIF ; cas general 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) z_s(ind) = interpol(z_zt(i_profil), s_z(i_profil), s_s(ind)) ; ; j'utilise la fonction spline pour interpoler le contenu ; l'interpolation lineaire marche aussi. Elle donne des resultats differents. Elle ; me semble moins propre. Je la donne en remarque. ; IF n_elements(i_profil) GT 2 THEN BEGIN c1_s(ind) = spline(z_zw(i_profil), c1_z(i_profil), z_s(ind), 1) ; c1_s(ind) = interpol(c1_z(i_profil), z_zw(i_profil), z_s(ind)) ENDIF ELSE BEGIN c1_s(ind) = interpol(c1_z(i_profil), z_zw(i_profil), z_s(ind)) ENDELSE ; ; x1_s(ind) = (c1_s(ind+1)-c1_s(ind))/(z_s(ind+1)-z_s(ind)) ENDIF ENDIF ; depth_bin (i, j, *) = z_s thick_bin (i, j, 0) = z_s(0) thick_bin (i, j, 1:N_s) = z_s(1:N_s)-z_s(0:N_s-1) x1_bin (i, j, *) = x1_s ; ENDFOR ENDFOR ; mask depth_bin with undefined values of x_bin depth_bin(where(finite(x1_bin, /nan) EQ 1)) = !VALUES.F_NAN thick_bin(where(finite(x1_bin, /nan) EQ 1)) = !VALUES.F_NAN END