source: trunk/SRC/ToBeReviewed/CALCULS/hdyn.pro @ 224

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
;
; @file_comments
; Calculate the height by rapport to a reference state for a depth reference.
; See keywords for different possibilities. By default, the state reference
; is rho=1020 and the depth reference is gdepw[ka] with ka the first W level
; directly above 1000 m.
;
; @categories
; Calculation
;
; @param TABSN {in}{required}
; array representing the salinity
;
; @param TABTN {in}{required}
; array representing the temperature. Has the same size than SN.
;
; @keyword GILL 
; We activate this key if we want to calculate the dynamic height
; like in the GILL page 215, which means by rapport to a reference state which
; vary in depth and which is determined by a reference temperature tref at 0°C
; and a reference salinity sref at 35 psu.
;
; @keyword LEVEL
; It is the same reference level to take. This level is defined like
; that gdepw[level] is the reference depth
;
; @keyword SREF
; Give a value to this keyword to change the reference salinity used in the
; calculation when GILL is activated.
; 
; @keyword TREF
; Give a value to this keyword to change the reference temperature used in the
; calculation when GILL is activated.
;
; @keyword PROFREF
; Give a depth to this keyword which will be considered as the reference depth
; (in this case, LEVEL has not any effect). the calculation will be effectuated
; until this depth effecting an interpolation between the the last W level above
; PROFREF and PROFREF.
;
; @keyword SURFACE_LEVEL {default=0}
; It is the level where we wan to calculate the dynamic height.
;
; @returns
; An array of the same size of sn and tn representing the dynamic height calculated
; from a reference depth nd by rapport to a reference state.
;
; @uses
; common.pro
;
; @restrictions
; Points for which we can not calculate the dynamic height (whose the batymetry
; is less deep than the reference depth) are put at the value !values.f_nan
;
; @restrictions
; approximation: The pressure in decibars is equal to the depth in meters 
; (the pressure increase of 1 bar every 10 m)
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION hdyn, tabsn, tabtn, TREF = tref, SREF = sref, PROFREF = profref, LEVEL = level, GILL = gill, SURFACE_LEVEL = surface_level
;
  compile_opt idl2, strictarrsubs
;
   tempsun = systime(1)         ; for key_performance
@common

   if NOT keyword_set(surface_level) then surface_level = 0
; useful if GILL is activated
   if NOT keyword_set(tref) then tref = 0.
   if NOT keyword_set(sref) then sref = 35.
; If needed, we determinate the reference depth and the W level situated directly above.
   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('arrays sn and tn must have the same size')
   if tailles[3] NE jpk then return, report('vertical dimension of sn and tn arrarrays must be equal to 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('not implemented')
   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