source: trunk/tools/density_binning/binning_neutral_and_co/bin_velocity_New.pro @ 161

FUNCTION bin_velocity_New, dens, sigma_values, x1, UU = uu, VV = vv, TT = tt
;
@common
;
; test the stability of dens and correct it if necessary
; ------------------------------------------------------
   dens=npc(dens)
;
; compute the density profile at U or V point if necessary
; --------------------------------------------------------   
   IF keyword_set(uu) THEN BEGIN
      density = boundperio(umask()*(dens+shift(dens, -1, 0, 0))/2, /uu, /orca2, /pscal)
      xmask = umask()
   ENDIF 
   IF keyword_set(vv) THEN BEGIN
      density = boundperio(vmask()*(dens+shift(dens, 0, -1, 0))/2, /vv, /orca2, /pscal)
      xmask = vmask()
   ENDIF 
      IF keyword_set(tt) THEN BEGIN
      density = dens
      xmask = tmask
   ENDIF 
;
; initializations
      s_s = sigma_values        ; density range over which we bin
      N_s = n_elements(s_s)     ; number of density levels
      N_z = jpk                 ; number of z levels
;   Profiles along z level ( suffixe _z)
      c1_z = fltarr(N_z)        ; profil du contenu vertical de x
      s_z = fltarr(N_z)         ; profil de la densite
; vertical coordinate: in z-coord only
      z_zt = gdept              ; profondeur au point T  (k=0 -> 5m)
                                ; z_zw = gdepw (old)       
      z_zw=shift(gdepw,-1)      ;                     W  (k=0 -> 10m)
      z_zw(jpk-1)=gdepw(jpk-1)  
      
; s or partial step case:
      ioMESH   = iodir+'meshes'
      fich_msh = 'mesh_mask_NDV1.nc'
      z3d_zt     = read_ncdf('e3v_ps',0,/timestep,iodir=ioMESH,/nostruct,/tout,filename=fich_msh)
      z3d_zt(where(z3d_zt eq 1.e20))=500.
      z3d_zw = total(z3d_zt,3,/cumulative)      ;  W  (k=0 -> 10m)
      z3d_zt = (shift(z3d_zw,0,0,1)+z3d_zw)*0.5
      z3d_zt(*,*,0) = z3d_zw(*,*,0)*0.5         ; T  (k=0 -> 5m)
      
;   Profiles along density level (suffixe _s)
      z_s  = fltarr(N_s+1)
      c1_s = fltarr(N_s+1)
      x1_s = fltarr(N_s+1)
;   x1 binned on density (output array)
      x1_bin = fltarr(jpi, jpj, N_s+1)
;
; vertical content of x1 (ensure  the integrale conservation)
      x1_content = fltarr(jpi, jpj, jpk)
      x1_content(*, *, 0) = e3t(0) * x1(*, *, 0) * xmask(*,*,0)
      FOR k = 1, jpk-1 DO x1_content(*, *, k) = x1_content(*, *, k-1) $ 
                                              + e3t(k) * x1(*, *, k) * xmask(*,*,k)
;
;  Loop over the horizontal domain 2D
;
   FOR i = 0, (jpi-1) DO BEGIN 
      FOR j = 0, (jpj-1) DO BEGIN 
      
; depth at T and W points
         z_zt(*) = z3d_zt(i,j,*)
         z_zw(*) = z3d_zw(i,j,*)
      
;
;  Indices des points T dans l ocean
;
         i_ocean = where(xmask(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
 
            s_z(*) = density(i, j, *)
            c1_z(*) = x1_content(i, j, *)
; critere supplementaire a imposer sur le profil pour eviter les cas
; pathologiques en attendant d'ecrire une vraie commande d'extraction
; de progils stritement croissant. Il s'agit donc d'un test adhoc.
            IF s_z(0) LT s_z(i_ocean(n_elements(i_ocean)-1)) THEN BEGIN 
;------------------------------------------------------------------------
; controle si le profil est bien strictement croissant (pour l'instant
; inutilisŽ (avis aux amateurs)
;
;            ds = (shift(s_z, -1)-s_z)(i_ocean(0:n_elements(i_ocean)-2))
;            ds(0) = 0
;            ind = where(ds LT 0.)
;            croissant =  ind(0) EQ -1
;------------------------------------------------------------------------
            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)
            
; 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
  ;gr          i_min = i_min(n_elements(i_min)-1)
  ;gr          i_max = i_max[0]
             i_min = i_min[0]
             i_max = i_max(n_elements(i_min)-1)
;            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 
;
     ds =(shift(s_z, -1)-s_z)(i_ocean(0:n_elements(i_ocean)-2))
     ind_c = where(ds LT 0.)     
;     croissant =  ind_c(0) EQ -1
;     IF ( i_min GT i_max ) then begin
;     print, 'bug ',i_min, i_max,'  en i,j=', i,j
;     endif
     IF ( ind_c(0) NE -1 ) then begin  
           print, 'bug   en i,j=', i,j,'  ind_c = ', ind_c
           stop
     endif
;     IF ( i_min LE i_max ) then begin
       IF ( ind_c(0) EQ -1 ) then begin
            ind = where( (s_s GE s_z(i_min)) AND (s_s LT s_z(i_max)) )
            IF ind[0] NE -1 THEN BEGIN 
               IF ( n_elements(ind) 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
;
                   c1_s(ind) = spline(z_zw(i_profil), c1_z(i_profil), z_s(ind), 1)
;
; 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))
;
; Remarque : on ne divise pas par l'epaisseur de la couche
               ENDIF
               ; cas pathologique (tout dans 1 gamme de densite)
               IF ( n_elements(ind) EQ 1 ) THEN BEGIN
 ;                  print, 'one density bin   en i,j=', i,j,'  ind = ', ind
                   c1_s(ind) = c1_z(jpk-1)
               ENDIF
;
                   x1_s(ind[0]-1) = (c1_s(ind[0])-c1_s(ind[0]-1));/(z_s(ind+1)-z_s(ind))
                   x1_s(ind) = (c1_s(ind+1)-c1_s(ind));/(z_s(ind+1)-z_s(ind))
;                   x1_s(ind[0]-1) = (c1_s(ind[0])-c1_s(ind[0]-1))/(z_s(ind+1)-z_s(ind))
;                   x1_s(ind) = (c1_s(ind+1)-c1_s(ind))/(z_s(ind+1)-z_s(ind))
            ENDIF 
     endif
         ENDIF 
         ENDIF
         x1_bin(i, j, *) = x1_s
         
      ENDFOR 
   ENDFOR 
   
   return, x1_bin
END