PRO bining3,  $                 ;
;INPUTS
   density,   $                 ; density at T points (3D array)
   x1,        $                 ; field at T points (e.g. Salinity or Temperature)
   sobwlmax,  $                 ; bowl depth array
   sig_bowl,  $                 ; switch for bowl overlay
;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)
   bowl_bin,  $                 ; bowl depth binned on density
    vol_bin2,  $                 ; volume of isopycnals 
; KEYWORDS
   SIGMA = sigma,      $        ; bining values
   DEPTH_T = depth_t,  $        ; depth of T level
   DEPTH_W = depth_w,  $        ; depth of W level
   E1T = e1t,          $
   E2T = e2t,          $
   E3T = e3t,          $
   TMASK = tmask                ; tmask (3D array)

   size3d = size(density)
   jpi = size3d(1)
   jpj = size3d(2)
   jpk = size3d(3)
   
   s_s = sigma
   N_s = n_elements(s_s)
   
   depth_bin  = fltarr(jpi, jpj, N_s)
   thick_bin  = fltarr(jpi, jpj, N_s)
   x1_bin = fltarr(jpi, jpj, N_s)
   bowl_bin = fltarr(jpi, jpj)
   vol_bin = fltarr(jpi, jpj, N_s)
   vol_bin2 = fltarr(N_s)

   x1_bin(*, *, *) = !VALUES.F_NAN 
   depth_bin(*, *, *) = !VALUES.F_NAN 
   thick_bin(*, *, *) = !VALUES.F_NAN 
   bowl_bin(*, *) = !VALUES.F_NAN 
   vol_bin(*, *, *) =0
   vol_bin2(*) =0

   delta_sig = s_s(1)-s_s(0)

   vol_tot = 0.
   count_bin = fltarr(jpi, jpj, N_s)
   count_bin2 = fltarr(N_s)
   count_bin(*, *, *) =0
   count_bin2(*) =0

; perform bining - loop over i,j,k

   FOR i = 0, (jpi-1) DO BEGIN 
      FOR j = 0, (jpj-1) DO BEGIN 
         i_ocean = where(tmask(i, j, *) EQ 1)
         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)
            FOR k = 0, i_bottom-1 DO BEGIN
               bin_index = floor((density(i, j, k)-s_s(0))/delta_sig+0.5)
               IF bin_index LT 0 OR bin_index GT (N_s-1) THEN BEGIN
;                  print, ' WARNING: bin index out of bounds - density(i, j, k), bin_index, N_s = ', density(i, j, k), bin_index, N_s
               ENDIF ELSE BEGIN 
                  vol_bin(i, j, bin_index) = vol_bin(i, j, bin_index) + e1t(i,j)*e2t(i,j)*e3t(k)
                  vol_bin2(bin_index) = vol_bin2(bin_index) + e1t(i,j)*e2t(i,j)*e3t(k)
                  vol_tot = vol_tot + e1t(i,j)*e2t(i,j)*e3t(k)
                  count_bin(i, j, bin_index) = count_bin(i, j, bin_index)+1
                  count_bin2(bin_index) = count_bin2(bin_index)+1
               ENDELSE  
            ENDFOR 
         ENDIF 
      ENDFOR 
   ENDFOR 
   
   mask_index = where(count_bin EQ 0)
   vol_bin(mask_index) = !VALUES.F_NAN
   mask_index2 = where(count_bin2 EQ 0)
;   vol_bin2(mask_index2) = !VALUES.F_NAN

   print, ' total volume of box (m3) = ', vol_tot
;   print, ' density bins = ', vol_bin2
;   print, ' binned volume (1e14 m3) = ', transpose(vol_bin)*1.e-14
   
END