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