[2] | 1 | ;------------------------------------------------------------ |
---|
| 2 | ;------------------------------------------------------------ |
---|
| 3 | ;------------------------------------------------------------ |
---|
| 4 | ;+ |
---|
| 5 | ; |
---|
[142] | 6 | ; @file_comments |
---|
[224] | 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. |
---|
[2] | 11 | ; |
---|
[142] | 12 | ; @categories |
---|
[157] | 13 | ; Calculation |
---|
[2] | 14 | ; |
---|
[142] | 15 | ; @param TABSN {in}{required} |
---|
| 16 | ; array representing the salinity |
---|
[2] | 17 | ; |
---|
[142] | 18 | ; @param TABTN {in}{required} |
---|
[224] | 19 | ; array representing the temperature. Has the same size than SN. |
---|
[2] | 20 | ; |
---|
[224] | 21 | ; @keyword GILL |
---|
| 22 | ; We activate this key if we want to calculate the dynamic height |
---|
| 23 | ; like in the GILL page 215, which means by rapport to a reference state which |
---|
| 24 | ; vary in depth and which is determined by a reference temperature tref at 0°C |
---|
| 25 | ; and a reference salinity sref at 35 psu. |
---|
[2] | 26 | ; |
---|
[142] | 27 | ; @keyword LEVEL |
---|
[224] | 28 | ; It is the same reference level to take. This level is defined like |
---|
[142] | 29 | ; that gdepw[level] is the reference depth |
---|
[224] | 30 | ; |
---|
| 31 | ; @keyword SREF |
---|
| 32 | ; Give a value to this keyword to change the reference salinity used in the |
---|
| 33 | ; calculation when GILL is activated. |
---|
[142] | 34 | ; |
---|
| 35 | ; @keyword TREF |
---|
[224] | 36 | ; Give a value to this keyword to change the reference temperature used in the |
---|
| 37 | ; calculation when GILL is activated. |
---|
| 38 | ; |
---|
[142] | 39 | ; @keyword PROFREF |
---|
[224] | 40 | ; Give a depth to this keyword which will be considered as the reference depth |
---|
| 41 | ; (in this case, LEVEL has not any effect). the calculation will be effectuated |
---|
| 42 | ; until this depth effecting an interpolation between the the last W level above |
---|
[142] | 43 | ; PROFREF and PROFREF. |
---|
[2] | 44 | ; |
---|
[163] | 45 | ; @keyword SURFACE_LEVEL {default=0} |
---|
| 46 | ; It is the level where we wan to calculate the dynamic height. |
---|
[2] | 47 | ; |
---|
[142] | 48 | ; @returns |
---|
[224] | 49 | ; An array of the same size of sn and tn representing the dynamic height calculated |
---|
[142] | 50 | ; from a reference depth nd by rapport to a reference state. |
---|
[2] | 51 | ; |
---|
[224] | 52 | ; @uses |
---|
[142] | 53 | ; common.pro |
---|
[2] | 54 | ; |
---|
[142] | 55 | ; @restrictions |
---|
[224] | 56 | ; Points for which we can not calculate the dynamic height (whose the batymetry |
---|
[142] | 57 | ; is less deep than the reference depth) are put at the value !values.f_nan |
---|
[2] | 58 | ; |
---|
[142] | 59 | ; @restrictions |
---|
[224] | 60 | ; approximation: The pressure in decibars is equal to the depth in meters |
---|
| 61 | ; (the pressure increase of 1 bar every 10 m) |
---|
[2] | 62 | ; |
---|
[142] | 63 | ; @history |
---|
[157] | 64 | ; Sebastien Masson (smasson\@lodyc.jussieu.fr) |
---|
[2] | 65 | ; |
---|
[142] | 66 | ; @version |
---|
| 67 | ; $Id$ |
---|
[2] | 68 | ; |
---|
| 69 | ;- |
---|
| 70 | ;------------------------------------------------------------ |
---|
| 71 | ;------------------------------------------------------------ |
---|
| 72 | ;------------------------------------------------------------ |
---|
[224] | 73 | FUNCTION hdyn, tabsn, tabtn, TREF = tref, SREF = sref, PROFREF = profref, LEVEL = level, GILL = gill, SURFACE_LEVEL = surface_level |
---|
[114] | 74 | ; |
---|
| 75 | compile_opt idl2, strictarrsubs |
---|
| 76 | ; |
---|
[142] | 77 | tempsun = systime(1) ; for key_performance |
---|
[2] | 78 | @common |
---|
| 79 | |
---|
| 80 | if NOT keyword_set(surface_level) then surface_level = 0 |
---|
[142] | 81 | ; useful if GILL is activated |
---|
[2] | 82 | if NOT keyword_set(tref) then tref = 0. |
---|
| 83 | if NOT keyword_set(sref) then sref = 35. |
---|
[142] | 84 | ; If needed, we determinate the reference depth and the W level situated directly above. |
---|
[2] | 85 | if keyword_set(profref) then begin |
---|
| 86 | rien = where(gdepw LE profref, level) |
---|
| 87 | level = level-1 |
---|
| 88 | za = gdepw[level] |
---|
[224] | 89 | ENDIF ELSE BEGIN |
---|
| 90 | if NOT keyword_set(level) then BEGIN |
---|
[2] | 91 | rien = where(gdepw LE 1000., level) |
---|
| 92 | level = level-1 |
---|
| 93 | ENDIF |
---|
| 94 | profref = gdepw[level] |
---|
| 95 | za = profref |
---|
[224] | 96 | ENDELSE |
---|
[2] | 97 | tailles = size(tabsn) |
---|
| 98 | taillet = size(tabtn) |
---|
| 99 | if total(tailles[0:tailles[0]] NE taillet[0:taillet[0]]) NE 0 then $ |
---|
[224] | 100 | return, report('arrays sn and tn must have the same size') |
---|
| 101 | if tailles[3] NE jpk then return, report('vertical dimension of sn and tn arrarrays must be equal to jpk') |
---|
[2] | 102 | nx = nxt |
---|
| 103 | ny = nyt |
---|
| 104 | case (size(tabsn))[0] OF |
---|
[224] | 105 | 3:BEGIN |
---|
[2] | 106 | case 1 of |
---|
| 107 | tailles[1] eq jpi and tailles[2] eq jpj: BEGIN |
---|
[25] | 108 | sn = tabsn[firstxt:lastxt, firstyt:lastyt, *] |
---|
| 109 | tn = tabtn[firstxt:lastxt, firstyt:lastyt, *] |
---|
[2] | 110 | end |
---|
[224] | 111 | tailles[1] eq nx and tailles[2] eq ny:BEGIN |
---|
[2] | 112 | sn = tabsn |
---|
| 113 | tn = tabtn |
---|
| 114 | end |
---|
| 115 | else:return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') |
---|
| 116 | ENDCASE |
---|
| 117 | if keyword_set(gill) then $ |
---|
| 118 | rhonref = rhon(replicate(sref,nx, ny, jpk),replicate(tref,nx, ny, jpk), /insitu) $ |
---|
| 119 | ELSE rhonref = 1020. |
---|
| 120 | vol=(rhonref-rhon(sn,tn, /insitu))/rhonref |
---|
| 121 | e33d = replicate(1, nx*ny)#e3t |
---|
| 122 | e33d = reform(e33d, nx, ny, jpk, /over) |
---|
[25] | 123 | terre = where(tmask[firstxt:lastxt, firstyt:lastyt, *] EQ 0) |
---|
[2] | 124 | if terre[0] NE -1 then vol[terre] = !values.f_nan |
---|
| 125 | case level of |
---|
| 126 | 0:hdyn =100.*(profref-gdepw[0])*vol[*, *, 0] |
---|
| 127 | 1:hdyn =100.*(vol*e33d)[*, *, 0]+(profref-gdepw[1])*vol[*, *, 1] |
---|
| 128 | ELSE:hdyn =100.*total( (vol*e33d)[*, *, surface_level: level-1], 3) $ |
---|
| 129 | +(profref-gdepw[level])*vol[*, *, level] |
---|
| 130 | endcase |
---|
| 131 | END |
---|
[224] | 132 | 4:BEGIN |
---|
[2] | 133 | case 1 of |
---|
| 134 | tailles[1] eq jpi and tailles[2] eq jpj AND tailles[4] EQ jpt: BEGIN |
---|
[25] | 135 | sn = tabsn[firstxt:lastxt, firstyt:lastyt, *, *] |
---|
| 136 | tn = tabtn[firstxt:lastxt, firstyt:lastyt, *, *] |
---|
[2] | 137 | end |
---|
[224] | 138 | tailles[1] eq nx and tailles[2] eq ny AND tailles[4] EQ jpt:BEGIN |
---|
[2] | 139 | sn = tabsn |
---|
| 140 | tn = tabtn |
---|
| 141 | end |
---|
| 142 | else:return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') |
---|
| 143 | endcase |
---|
| 144 | if keyword_set(gill) then $ |
---|
| 145 | rhonref = rhon(replicate(sref,nx, ny, jpk, jpt),replicate(tref,nx, ny, jpk, jpt), /insitu) $ |
---|
| 146 | ELSE rhonref = 1020. |
---|
| 147 | vol=(rhonref-rhon(sn,tn, /insitu))/rhonref |
---|
| 148 | e33d = replicate(1, nx*ny)#e3t |
---|
| 149 | e33d = e33d[*]#replicate(1, jpt) |
---|
| 150 | e33d = reform(e33d, nx, ny, jpk, jpt, /over) |
---|
[25] | 151 | mask = tmask[firstxt:lastxt, firstyt:lastyt, *] |
---|
[2] | 152 | mask = mask[*]#replicate(1, jpt) |
---|
| 153 | terre = where(mask EQ 0) |
---|
| 154 | if terre[0] NE -1 then vol[terre] = !values.f_nan |
---|
| 155 | case level of |
---|
| 156 | 0:hdyn =100.*(profref-gdepw[0])*vol[*, *, 0, *] |
---|
| 157 | 1:hdyn =100.*(vol*e33d)[*, *, 0, *]+(profref-gdepw[1])*vol[*, *, 1, *] |
---|
| 158 | ELSE:hdyn =100.*total( (vol*e33d)[*, *, surface_level: level-1, *], 3) $ |
---|
| 159 | +(profref-gdepw[level])*vol[*, *, level, *] |
---|
| 160 | endcase |
---|
| 161 | END |
---|
[224] | 162 | ELSE: return, report('not implemented') |
---|
[2] | 163 | ENDCASE |
---|
| 164 | varunit = 'cm' |
---|
| 165 | varname = 'Dynamic Height (href='+strtrim(round(profref), 1)+'m)' |
---|
[224] | 166 | IF keyword_set(key_performance) THEN print, 'temps hdyn', systime(1)-tempsun |
---|
[2] | 167 | |
---|
| 168 | return, hdyn |
---|
| 169 | end |
---|