source: trunk/SRC/ToBeReviewed/CALCULS/projectondepth.pro

;+
;
; @file_comments
; Allows to project a 3d field following a depth array.
;
; @categories
; Without loop
;
; @param ARRAYIN {type=3d array}
; It is a 3d array whose 3rd dimension must be equal to jpk
;
; @param DEPTHIN {type=2d array}
; It is a 2d array indicating for each point n, at which depth to project
;
; @returns
; A 2d array which is the projection of the 3d array following depths indicated by depthin
;
; @uses
; <pro>common</pro>
;
; @restrictions
; points at !values.f_nan impossible calculation. Land points masked at valmask.
;
; @examples
; we build a possible depths array
;   IDL> a=gdept[jpk-1]/(1.*jpi*jpj)*findgen(jpi,jpj)
; We build an array to project on these depths. For the test,
; we build a 3d array whose each vector following z is the depth.
;   IDL> arraytest=replicate(1,jpi*jpj)#gdept
;   IDL> arraytest=reform(arraytest,jpi,jpj,jpk, /over)
; We test the projection of the depth array on the depth...
;   IDL> plt, 1e6*(a-projectondepth(arraytest,a)),/nocontour
;   ->null field at 1e-6 pres
;
;  verification projecting the temperature of 20°C for example...
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                      15/6/2000
;
; @version
; $Id$
;
;-
FUNCTION projectondepth, arrayin, depthin
;
  compile_opt idl2, strictarrsubs
;
   tempsun = systime(1)         ; To key_performance
@common
;------------------------------------------------------------
   depth = litchamp(depthin)
   array = litchamp(arrayin)
; Small 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 of the coherence between array's size and the domain
   grille, mask, -1, -1, -1,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz
   case 1 of
      tailledepth[1] eq jpi and tailledepth[2] eq jpj:depth=depth[firstx:lastx, firsty:lasty]
      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('2d array must have its 3d dimension equal to jpk')
      taillearray[1] eq jpi and taillearray[2] eq jpj:array=array[firstx:lastx, firsty:lasty, *]
      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)
; we delete points at !values.f_nan
   notanumber = where(finite(flevel, /nan) EQ 1)
   if notanumber[0] NE -1 then flevel[notanumber] = 0
; we sill (delete land points at valmask for example)
   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)
;
; We put back points at !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
; We mask land points at valmask
   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