;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME:projectondepth
;
; PURPOSE: routine permettant de projeter un champ 3d suivant un
; tableau de profondeurs.
;
; CATEGORY: sans boucles
;
; CALLING SEQUENCE:res=projectondepth(arrayin, depthin)
; 
; INPUTS:
;       arrayin: un tableau 3d dont la 3eme dimension doit etre egale
;       a jpk
;       depthin: un tableau 2d indiquant n chaque point a quel
;       profondeur projeter
;       
; KEYWORD PARAMETERS:none
;
; OUTPUTS:res: un tableau 2d projection du tableau 3d suivant les
; profondeurs indiquees par depthin
;
; COMMON BLOCKS:common.pro
;
; SIDE EFFECTS: points a !values.f_nan qd calcul impossible. points
; terres masques a Valmask.
;
; RESTRICTIONS:
;
; EXAMPLE:
;
; on contruit un tableau de profondeurs possibles
;   IDL> a=gdept[jpk-1]/(1.*jpi*jpj)*findgen(jpi,jpj)
; on contruit un tableau a projeter sur ces profondeurs. pour le test
; on construit un tableau 3d dont chaque vecteur suivant z est la
; profondeur. 
;   IDL> arraytest=replicate(1,jpi*jpj)#gdept
;   IDL> arraytest=reform(arraytest,jpi,jpj,jpk, /over)
; on test la projection du tabeau profondeur sur la profondeur...
;   IDL> plt, 1e6*(a-projectondepth(arraytest,a)),/nocontour
;   ->champ nul a 1e-6 pres
; 
;  verifcation en projettant la temperature sur la profondeur de la 2o
;  degres par exemple...
;
; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
;                      15/6/2000
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION projectondepth, arrayin, depthin
   tempsun = systime(1)         ; pour key_performance
@common
;------------------------------------------------------------
   depth = litchamp(depthin)
   array = litchamp(arrayin)
; petites verifications
   tailledepth = size(depth)
   taillearray = size(array)
   if tailledepth[0] NE 2 THEN return, report('Depth array must have 2 dimensions')
   if taillearray[0] NE 3 THEN return, report('Array in must have 3 dimensions')
; verification de la coherence entre la taille du tableau et le
; domaine 
   grille,mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz
   case 1 of
      tailledepth[1] eq jpi and tailledepth[2] eq jpj:depth=depth[premierx:dernierx, premiery:derniery]
      tailledepth[1] eq  nx and tailledepth[2] eq  ny:
      else:return, report('Probleme d''adequation entre les tailles du domaine et celle du tableau de profondeur')
   endcase
   case 1 OF
      taillearray[3] NE jpk:return, report('Le tableau 3d doit avoir sa 3eme dimension egale a jpk')
      taillearray[1] eq jpi and taillearray[2] eq jpj:array=array[premierx:dernierx, premiery:derniery, *]
      taillearray[1] eq  nx and taillearray[2] eq  ny:
      else:return, report('Probleme d''adequation entre les tailles du domaine et celle du tableau de profondeur')
   endcase
;
; c''est parti
;
   flevel = depth2floatlevel(depth)
; on vire les points a !values.f_nan
   notanumber = where(finite(flevel, /nan) EQ 1)
   if notanumber[0] NE -1 then flevel[notanumber] = 0
; on seuil (vire les points terres a valmask par ex)
   flevel = 0 > flevel < (jpk-1)
;
   indexup = level2index(floor(flevel))
   indexlow = nx*ny+indexup
   out = where(indexlow GE nx*ny*jpk-1)
   if out[0] NE -1 then indexlow[out] = indexlow[out]-nx*ny
;
   weight = flevel-floor(flevel)
   res = array[indexup]
   res = res+weight*(array[indexlow]-res)
;
; on replace les points a !values.f_nan
   if notanumber[0] NE -1 then res[notanumber] = !values.f_nan
   if out[0] NE -1 then res[out] = !values.f_nan
; on masque les points terres a valmask
   grille,mask
   if n_elements(valmask) EQ 0 then valmask = 1e20
   terre = where((temporary(mask))[*, *, 0] EQ 0)
   if terre[0] NE -1 then res[terre] = valmask
;------------------------------------------------------------
   if keyword_set(key_performance) THEN print, 'temps projectondepth', systime(1)-tempsun 
   return, res
end