source: trunk/ToBeReviewed/CALCULS/hdyn.pro @ 25

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME:hdyn
;
; PURPOSE:calcule la hauteur dynamique par rapport a un etat de
; reference pour une profondeur de reference. Cf les mots cles pour
; les differentes possibilites. 
; Par defaut l''etat de reference est
; rho=1020 et la profondeur de reference est gdepw[ka] avec ka le
; premier niveau W directement au dessus de 1000 m. 
;
; CATEGORY: calculs de post-traitement
;
; CALLING SEQUENCE:res=hdyn(sn,tn)
;
; INPUTS:sn et tn sont des tableaux de meme taille representant la
; salinite et la temperature.
;
; KEYWORD PARAMETERS:
;
;        GILL: activer cette cle si on veut faire le calcul de la
;        hauteur dynamique comme ds le GILL page 215 cad par rapport a
;        un etat de reference qui varie en profondeur et qui est
;        determine par une temperature de reference tref a 0 degre et
;        une salinite de reference sref a 35psu
;
;        LEVEL: C''est le niveau de reference a prendre. Ce niveau est
;        definit tel que gdepw[level] est la profondeur de reference
;
;        SREF: donner une valeur a ce mot cle pour changer la salinite
;        de reference utiliser ds le calcul lorsque GILL est active.
;   
;        TREF: donner une valeur a ce mot cle pour changer la temperature
;        de reference utiliser ds le calcul lorsque GILL est active.
;
;        PROFREF: donner a ce mot cle une profondeur qui sera prise comme
;        la profondeur de reference (ds ce cas LEVEL n'' a aucun
;        effet). le calcul sera alors effectue jusqu''a cette
;        profondeur en effectuant une interpolation entre le dernier
;        niveau W au dessus de PROFREF et PROFREF.
;
;        SURFACE_LEVEL: C''est le niveau auquel on veut calculer la
;        hauteur dynamique. Par defaut c''est le niveau 0.
;
; OUTPUTS:un tableau de la meme taille que sn et tn representant la
; hauteur dynamique calculee a partir d''une profondeur de reference
; et par rapport a un etat de reference.
;
; COMMON BLOCKS:
;       common.pro
;
; SIDE EFFECTS:
; les points pour lesquels on nje peut calcule la hauteur dynamique
; (dont la batymetrie est moins profonde que la profondeur de
; reference) sont mis a la valeur !values.f_nan
;
; RESTRICTIONS:approximation: la pression en decibars est egale a la
; profondeur en metres (la pression augmente de 1bar tous les 10m)
;
; EXAMPLE:
;
; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
;
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION hdyn,  tabsn, tabtn, TREF = tref,  SREF = sref, PROFREF = profref, LEVEL = level, GILL = gill, SURFACE_LEVEL = surface_level
   tempsun = systime(1)         ; pour key_performance
@common

   if NOT keyword_set(surface_level) then surface_level = 0
; utile si GILL est active
   if NOT keyword_set(tref) then tref = 0.
   if NOT keyword_set(sref) then sref = 35.
; on determine si besoin est la profondeur de reference et le niveau
; W situe directement au dessus.
   if keyword_set(profref) then begin
      rien = where(gdepw LE profref, level)
      level = level-1
      za = gdepw[level]
   ENDIF ELSE BEGIN 
      if NOT keyword_set(level) then BEGIN 
         rien = where(gdepw LE 1000., level)
         level = level-1
      ENDIF
      profref = gdepw[level]
      za = profref
   ENDELSE 
   tailles = size(tabsn)
   taillet = size(tabtn)
   if total(tailles[0:tailles[0]] NE taillet[0:taillet[0]]) NE 0 then $
    return,  report('Les tableaux sn et tn doivent avoir la meme taille')
   if tailles[3] NE jpk then return, report('La dim verticale des tableaux sn et tn doit etre egalre a jpk')
   nx = nxt
   ny = nyt
   case (size(tabsn))[0] OF
      3:BEGIN 
         case 1 of
            tailles[1] eq jpi and tailles[2] eq jpj: BEGIN
               sn = tabsn[firstxt:lastxt, firstyt:lastyt, *]
               tn = tabtn[firstxt:lastxt, firstyt:lastyt, *]
            end
            tailles[1] eq  nx and tailles[2] eq  ny:BEGIN 
               sn = tabsn
               tn = tabtn
            end
            else:return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
         ENDCASE
         if keyword_set(gill) then $
          rhonref = rhon(replicate(sref,nx, ny, jpk),replicate(tref,nx, ny, jpk), /insitu) $
         ELSE rhonref = 1020.
         vol=(rhonref-rhon(sn,tn, /insitu))/rhonref
         e33d = replicate(1, nx*ny)#e3t
         e33d = reform(e33d, nx, ny, jpk, /over)
         terre = where(tmask[firstxt:lastxt, firstyt:lastyt, *] EQ 0)
         if terre[0] NE -1 then vol[terre] = !values.f_nan
         case level of
            0:hdyn =100.*(profref-gdepw[0])*vol[*, *, 0]
            1:hdyn =100.*(vol*e33d)[*, *, 0]+(profref-gdepw[1])*vol[*, *, 1]
            ELSE:hdyn =100.*total( (vol*e33d)[*, *, surface_level: level-1], 3) $
             +(profref-gdepw[level])*vol[*, *, level]
         endcase
      END
      4:BEGIN 
         case 1 of
            tailles[1] eq jpi and tailles[2] eq jpj AND tailles[4] EQ jpt: BEGIN
               sn = tabsn[firstxt:lastxt, firstyt:lastyt, *, *]
               tn = tabtn[firstxt:lastxt, firstyt:lastyt, *, *]
            end
            tailles[1] eq  nx and tailles[2] eq  ny AND tailles[4] EQ jpt:BEGIN 
               sn = tabsn
               tn = tabtn
            end
            else:return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
         endcase
         if keyword_set(gill) then $
          rhonref = rhon(replicate(sref,nx, ny, jpk, jpt),replicate(tref,nx, ny, jpk, jpt), /insitu) $
         ELSE rhonref = 1020.
         vol=(rhonref-rhon(sn,tn, /insitu))/rhonref
         e33d = replicate(1, nx*ny)#e3t
         e33d = e33d[*]#replicate(1, jpt)
         e33d = reform(e33d, nx, ny, jpk, jpt, /over)
         mask = tmask[firstxt:lastxt, firstyt:lastyt, *]
         mask = mask[*]#replicate(1, jpt)
         terre = where(mask EQ 0)
         if terre[0] NE -1 then vol[terre] = !values.f_nan
         case level of
            0:hdyn =100.*(profref-gdepw[0])*vol[*, *, 0, *]
            1:hdyn =100.*(vol*e33d)[*, *, 0, *]+(profref-gdepw[1])*vol[*, *, 1, *]
            ELSE:hdyn =100.*total( (vol*e33d)[*, *, surface_level: level-1, *], 3) $
             +(profref-gdepw[level])*vol[*, *, level, *]
         endcase
      END
      ELSE: return,  report('cas non code')
   ENDCASE
   varunit = 'cm'
   varname = 'Dynamic Height (href='+strtrim(round(profref), 1)+'m)'
   IF keyword_set(key_performance) THEN print, 'temps hdyn', systime(1)-tempsun 

   return, hdyn
end