1 | ;------------------------------------------------------------ |
---|
2 | ;------------------------------------------------------------ |
---|
3 | ;------------------------------------------------------------ |
---|
4 | ;+ |
---|
5 | ; |
---|
6 | ; @file_comments |
---|
7 | ; Calculate the height by rapport to a reference state for a depth reference. |
---|
8 | ; See keywords for different possibilities. By default, the state reference |
---|
9 | ; is rho=1020 and the depth reference is gdepw[ka] with ka the first W level |
---|
10 | ; directly above 1000 m. |
---|
11 | ; |
---|
12 | ; @categories |
---|
13 | ; Calculation |
---|
14 | ; |
---|
15 | ; @param TABSN {in}{required} |
---|
16 | ; array representing the salinity |
---|
17 | ; |
---|
18 | ; @param TABTN {in}{required} |
---|
19 | ; array representing the temperature.Has the same size than SN |
---|
20 | ; |
---|
21 | ; @keyword GILL We activate this key if we want to calculate the dynamic height |
---|
22 | ; like in the GILL page 215, which means by rapport to a reference state which |
---|
23 | ; vary in depth and which is determined by a reference temperature tref at 0°C |
---|
24 | ; and a reference salinity sref at 35psu |
---|
25 | ; |
---|
26 | ; @keyword LEVEL |
---|
27 | ; It is the same reference level to take. This level is defined like |
---|
28 | ; that gdepw[level] is the reference depth |
---|
29 | ; |
---|
30 | ; @keyword SREF |
---|
31 | ; Give a value to this keyword to change the reference salinity used in the |
---|
32 | ; calculation when GILL is activated. |
---|
33 | ; |
---|
34 | ; @keyword TREF |
---|
35 | ; Give a value to this keyword to change the reference temperature used in the |
---|
36 | ; calculation when GILL is activated. |
---|
37 | ; |
---|
38 | ; @keyword PROFREF |
---|
39 | ; Give a depth to this keyword which will be considered as the reference depth |
---|
40 | ; (in this case, LEVEL has not any effect). the calculation will be effectuated |
---|
41 | ; until this depth effecting an interpolation between the the last W level above |
---|
42 | ; PROFREF and PROFREF. |
---|
43 | ; |
---|
44 | ; @keyword SURFACE_LEVEL {default=0} |
---|
45 | ; It is the level where we wan to calculate the dynamic height. |
---|
46 | ; |
---|
47 | ; @returns |
---|
48 | ; An array of the same size of sn and tn representing the dynamic height calculated |
---|
49 | ; from a reference depth nd by rapport to a reference state. |
---|
50 | ; |
---|
51 | ; @uses |
---|
52 | ; common.pro |
---|
53 | ; |
---|
54 | ; @restrictions |
---|
55 | ; Points for which we can not calculate the dynamic height (whose the batymetry |
---|
56 | ; is less deep than the reference depth) are put at the value !values.f_nan |
---|
57 | ; |
---|
58 | ; @restrictions |
---|
59 | ; approximation: The pressure in decibars is equal to the depth in meters (the pressure increase of 1bar all 10m) |
---|
60 | ; |
---|
61 | ; @history |
---|
62 | ; Sebastien Masson (smasson\@lodyc.jussieu.fr) |
---|
63 | ; |
---|
64 | ; @version |
---|
65 | ; $Id$ |
---|
66 | ; |
---|
67 | ;- |
---|
68 | ;------------------------------------------------------------ |
---|
69 | ;------------------------------------------------------------ |
---|
70 | ;------------------------------------------------------------ |
---|
71 | FUNCTION hdyn, tabsn, tabtn, TREF = tref, SREF = sref, PROFREF = profref, LEVEL = level, GILL = gill, SURFACE_LEVEL = surface_level |
---|
72 | ; |
---|
73 | compile_opt idl2, strictarrsubs |
---|
74 | ; |
---|
75 | tempsun = systime(1) ; for key_performance |
---|
76 | @common |
---|
77 | |
---|
78 | if NOT keyword_set(surface_level) then surface_level = 0 |
---|
79 | ; useful if GILL is activated |
---|
80 | if NOT keyword_set(tref) then tref = 0. |
---|
81 | if NOT keyword_set(sref) then sref = 35. |
---|
82 | ; If needed, we determinate the reference depth and the W level situated directly above. |
---|
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 |
---|
167 | end |
---|