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

Last change on this file since 114 was 114, checked in by smasson, 18 years ago

new compilation options (compile_opt idl2, strictarrsubs) in each routine

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 6.7 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;+
5; NAME:hdyn
6;
7; PURPOSE:calcule la hauteur dynamique par rapport a un etat de
8; reference pour une profondeur de reference. Cf les mots cles pour
9; les differentes possibilites.
10; Par defaut l''etat de reference est
11; rho=1020 et la profondeur de reference est gdepw[ka] avec ka le
12; premier niveau W directement au dessus de 1000 m.
13;
14; CATEGORY: calculs de post-traitement
15;
16; CALLING SEQUENCE:res=hdyn(sn,tn)
17;
18; INPUTS:sn et tn sont des tableaux de meme taille representant la
19; salinite et la temperature.
20;
21; KEYWORD PARAMETERS:
22;
23;        GILL: activer cette cle si on veut faire le calcul de la
24;        hauteur dynamique comme ds le GILL page 215 cad par rapport a
25;        un etat de reference qui varie en profondeur et qui est
26;        determine par une temperature de reference tref a 0 degre et
27;        une salinite de reference sref a 35psu
28;
29;        LEVEL: C''est le niveau de reference a prendre. Ce niveau est
30;        definit tel que gdepw[level] est la profondeur de reference
31;
32;        SREF: donner une valeur a ce mot cle pour changer la salinite
33;        de reference utiliser ds le calcul lorsque GILL est active.
34;   
35;        TREF: donner une valeur a ce mot cle pour changer la temperature
36;        de reference utiliser ds le calcul lorsque GILL est active.
37;
38;        PROFREF: donner a ce mot cle une profondeur qui sera prise comme
39;        la profondeur de reference (ds ce cas LEVEL n'' a aucun
40;        effet). le calcul sera alors effectue jusqu''a cette
41;        profondeur en effectuant une interpolation entre le dernier
42;        niveau W au dessus de PROFREF et PROFREF.
43;
44;        SURFACE_LEVEL: C''est le niveau auquel on veut calculer la
45;        hauteur dynamique. Par defaut c''est le niveau 0.
46;
47; OUTPUTS:un tableau de la meme taille que sn et tn representant la
48; hauteur dynamique calculee a partir d''une profondeur de reference
49; et par rapport a un etat de reference.
50;
51; COMMON BLOCKS:
52;       common.pro
53;
54; SIDE EFFECTS:
55; les points pour lesquels on nje peut calcule la hauteur dynamique
56; (dont la batymetrie est moins profonde que la profondeur de
57; reference) sont mis a la valeur !values.f_nan
58;
59; RESTRICTIONS:approximation: la pression en decibars est egale a la
60; profondeur en metres (la pression augmente de 1bar tous les 10m)
61;
62; EXAMPLE:
63;
64; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
65;
66;-
67;------------------------------------------------------------
68;------------------------------------------------------------
69;------------------------------------------------------------
70FUNCTION hdyn,  tabsn, tabtn, TREF = tref,  SREF = sref, PROFREF = profref, LEVEL = level, GILL = gill, SURFACE_LEVEL = surface_level
71;
72  compile_opt idl2, strictarrsubs
73;
74   tempsun = systime(1)         ; pour key_performance
75@common
76
77   if NOT keyword_set(surface_level) then surface_level = 0
78; utile si GILL est active
79   if NOT keyword_set(tref) then tref = 0.
80   if NOT keyword_set(sref) then sref = 35.
81; on determine si besoin est la profondeur de reference et le niveau
82; W situe directement au dessus.
83   if keyword_set(profref) then begin
84      rien = where(gdepw LE profref, level)
85      level = level-1
86      za = gdepw[level]
87   ENDIF ELSE BEGIN
88      if NOT keyword_set(level) then BEGIN
89         rien = where(gdepw LE 1000., level)
90         level = level-1
91      ENDIF
92      profref = gdepw[level]
93      za = profref
94   ENDELSE
95   tailles = size(tabsn)
96   taillet = size(tabtn)
97   if total(tailles[0:tailles[0]] NE taillet[0:taillet[0]]) NE 0 then $
98    return,  report('Les tableaux sn et tn doivent avoir la meme taille')
99   if tailles[3] NE jpk then return, report('La dim verticale des tableaux sn et tn doit etre egalre a jpk')
100   nx = nxt
101   ny = nyt
102   case (size(tabsn))[0] OF
103      3:BEGIN
104         case 1 of
105            tailles[1] eq jpi and tailles[2] eq jpj: BEGIN
106               sn = tabsn[firstxt:lastxt, firstyt:lastyt, *]
107               tn = tabtn[firstxt:lastxt, firstyt:lastyt, *]
108            end
109            tailles[1] eq  nx and tailles[2] eq  ny:BEGIN
110               sn = tabsn
111               tn = tabtn
112            end
113            else:return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
114         ENDCASE
115         if keyword_set(gill) then $
116          rhonref = rhon(replicate(sref,nx, ny, jpk),replicate(tref,nx, ny, jpk), /insitu) $
117         ELSE rhonref = 1020.
118         vol=(rhonref-rhon(sn,tn, /insitu))/rhonref
119         e33d = replicate(1, nx*ny)#e3t
120         e33d = reform(e33d, nx, ny, jpk, /over)
121         terre = where(tmask[firstxt:lastxt, firstyt:lastyt, *] EQ 0)
122         if terre[0] NE -1 then vol[terre] = !values.f_nan
123         case level of
124            0:hdyn =100.*(profref-gdepw[0])*vol[*, *, 0]
125            1:hdyn =100.*(vol*e33d)[*, *, 0]+(profref-gdepw[1])*vol[*, *, 1]
126            ELSE:hdyn =100.*total( (vol*e33d)[*, *, surface_level: level-1], 3) $
127             +(profref-gdepw[level])*vol[*, *, level]
128         endcase
129      END
130      4:BEGIN
131         case 1 of
132            tailles[1] eq jpi and tailles[2] eq jpj AND tailles[4] EQ jpt: BEGIN
133               sn = tabsn[firstxt:lastxt, firstyt:lastyt, *, *]
134               tn = tabtn[firstxt:lastxt, firstyt:lastyt, *, *]
135            end
136            tailles[1] eq  nx and tailles[2] eq  ny AND tailles[4] EQ jpt:BEGIN
137               sn = tabsn
138               tn = tabtn
139            end
140            else:return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
141         endcase
142         if keyword_set(gill) then $
143          rhonref = rhon(replicate(sref,nx, ny, jpk, jpt),replicate(tref,nx, ny, jpk, jpt), /insitu) $
144         ELSE rhonref = 1020.
145         vol=(rhonref-rhon(sn,tn, /insitu))/rhonref
146         e33d = replicate(1, nx*ny)#e3t
147         e33d = e33d[*]#replicate(1, jpt)
148         e33d = reform(e33d, nx, ny, jpk, jpt, /over)
149         mask = tmask[firstxt:lastxt, firstyt:lastyt, *]
150         mask = mask[*]#replicate(1, jpt)
151         terre = where(mask EQ 0)
152         if terre[0] NE -1 then vol[terre] = !values.f_nan
153         case level of
154            0:hdyn =100.*(profref-gdepw[0])*vol[*, *, 0, *]
155            1:hdyn =100.*(vol*e33d)[*, *, 0, *]+(profref-gdepw[1])*vol[*, *, 1, *]
156            ELSE:hdyn =100.*total( (vol*e33d)[*, *, surface_level: level-1, *], 3) $
157             +(profref-gdepw[level])*vol[*, *, level, *]
158         endcase
159      END
160      ELSE: return,  report('cas non code')
161   ENDCASE
162   varunit = 'cm'
163   varname = 'Dynamic Height (href='+strtrim(round(profref), 1)+'m)'
164   IF keyword_set(key_performance) THEN print, 'temps hdyn', systime(1)-tempsun
165
166   return, hdyn
167end
Note: See TracBrowser for help on using the repository browser.