source: trunk/procs/bin @ 13

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