| 1 | ;+ |
|---|
| 2 | ; |
|---|
| 3 | ; @file_comments |
|---|
| 4 | ; calculate the norm of vectors field located on Arakawa C-grid |
|---|
| 5 | ; |
|---|
| 6 | ; @categories |
|---|
| 7 | ; Calculation |
|---|
| 8 | ; |
|---|
| 9 | ; @param UU {in}{required} |
|---|
| 10 | ; Matrix representing the zonal coordinates (at U/V point) of a field of vectors |
|---|
| 11 | ; A 2D (xy), 3D (xyz or yt), 4D (xyzt) or a structure readable by |
|---|
| 12 | ; <pro>litchamp</pro> and containing a 2D (xy), 3D (xyz or yt), 4D (xyzt) array. |
|---|
| 13 | ; Note that the dimension of the array must suit the domain dimension. |
|---|
| 14 | ; |
|---|
| 15 | ; @param VV {in}{required} |
|---|
| 16 | ; Matrix representing the meridional coordinates (at V/U point) of a field of vectors |
|---|
| 17 | ; A 2D (xy), 3D (xyz or yt), 4D (xyzt) or a structure readable by |
|---|
| 18 | ; <pro>litchamp</pro> and containing a 2D (xy), 3D (xyz or yt), 4D (xyzt) array. |
|---|
| 19 | ; Note that the dimension of the array must suit the domain dimension. |
|---|
| 20 | ; |
|---|
| 21 | ; @keyword DIREC |
|---|
| 22 | ; 't' 'x' 'y' 'z' 'xys' 'xz' 'yz' 'xyz' 'xt' 'yt' 'zt' 'xyt' |
|---|
| 23 | ; 'xzt' 'yzt' 'xyzt' Direction on which do averages |
|---|
| 24 | ; |
|---|
| 25 | ; @keyword _EXTRA |
|---|
| 26 | ; Used to declare that this routine accepts inherited keywords |
|---|
| 27 | ; |
|---|
| 28 | ; @returns |
|---|
| 29 | ; A 2D (xy), 3D (xyz or yt), 4D (xyzt) Array |
|---|
| 30 | ; |
|---|
| 31 | ; @uses |
|---|
| 32 | ; @cm_4mesh |
|---|
| 33 | ; @cm_4data |
|---|
| 34 | ; @cm_4cal |
|---|
| 35 | ; |
|---|
| 36 | ; @restrictions |
|---|
| 37 | ; The norm is calculated on points T. To do this calculation, we average |
|---|
| 38 | ; field U and V on points T before calculate the norm. At the edge of |
|---|
| 39 | ; coast and of domain, we can not calculate fields U and V at points T, |
|---|
| 40 | ; that is why these points are at value !values.f_nan. |
|---|
| 41 | ; |
|---|
| 42 | ; When we calculate on a reduce geographic domain, field U and V have not |
|---|
| 43 | ; necessarily the same number of point. In this case, we recut U and V to |
|---|
| 44 | ; keep only common points. We profit of this to redo a <pro>domdef</pro> which redefine |
|---|
| 45 | ; a geographic domain on which fields U and V are extracted on same points |
|---|
| 46 | ; |
|---|
| 47 | ; To know what type of array we work with, we test its size and dates |
|---|
| 48 | ; gave by time[0] and time[jpt-1] to know if thee is a temporal dimension. |
|---|
| 49 | ; Before to start norm, make sure that time and jpt are defined how |
|---|
| 50 | ; they have to! |
|---|
| 51 | ; |
|---|
| 52 | ; @examples |
|---|
| 53 | ; To calculate the average of the norm of streams on all the domain |
|---|
| 54 | ; between 0 and 50: |
|---|
| 55 | ; IDL> domdef, 0, 50 |
|---|
| 56 | ; IDL> res = norm(un, vn, dir = 'xyz') |
|---|
| 57 | ; |
|---|
| 58 | ; @history |
|---|
| 59 | ; Sebastien Masson (smasson\@lodyc.jussieu.fr) |
|---|
| 60 | ; 9/6/1999 |
|---|
| 61 | ; @version |
|---|
| 62 | ; $Id$ |
|---|
| 63 | ; |
|---|
| 64 | ;- |
|---|
| 65 | ; |
|---|
| 66 | FUNCTION norm, uu, vv, DIREC = direc, _EXTRA = ex |
|---|
| 67 | ; |
|---|
| 68 | compile_opt idl2, strictarrsubs |
|---|
| 69 | ; |
|---|
| 70 | @cm_4mesh |
|---|
| 71 | @cm_4data |
|---|
| 72 | @cm_4cal |
|---|
| 73 | IF NOT keyword_set(key_forgetold) THEN BEGIN |
|---|
| 74 | @updatenew |
|---|
| 75 | @updatekwd |
|---|
| 76 | ENDIF |
|---|
| 77 | ;--------------------------------------------------------- |
|---|
| 78 | time1 = systime(1) ; To key_performance |
|---|
| 79 | ; |
|---|
| 80 | IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ |
|---|
| 81 | return, report(['This version of norm is based on Arakawa C-grid.' $ |
|---|
| 82 | , 'U and V grids must therefore be defined']) |
|---|
| 83 | ;------------------------------------------------------------ |
|---|
| 84 | if NOT keyword_set(direc) then direc = 0 |
|---|
| 85 | ; construction of u and v at points T |
|---|
| 86 | u = litchamp(uu) |
|---|
| 87 | v = litchamp(vv) |
|---|
| 88 | date1 = time[0] |
|---|
| 89 | if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1] |
|---|
| 90 | |
|---|
| 91 | if (size(u))[0] NE (size(v))[0] then return, -1 |
|---|
| 92 | ; |
|---|
| 93 | grilleu = litchamp(uu, /grid) |
|---|
| 94 | if grilleu EQ '' then grilleu = 'U' |
|---|
| 95 | grillev = litchamp(vv, /grid) |
|---|
| 96 | if grillev EQ '' then grillev = 'V' |
|---|
| 97 | IF grilleu EQ 'V' AND grillev EQ 'U' THEN inverse = 1 |
|---|
| 98 | IF grilleu EQ 'T' AND grillev EQ 'T' THEN BEGIN |
|---|
| 99 | interpolle = 0 |
|---|
| 100 | return, report('Case not coded, but easy to do...!') |
|---|
| 101 | ENDIF ELSE interpolle = 1 |
|---|
| 102 | if keyword_set(inverse) then begin |
|---|
| 103 | tmp = u |
|---|
| 104 | u = temporary(v) |
|---|
| 105 | v = temporary(tmp) |
|---|
| 106 | endif |
|---|
| 107 | ;------------------------------------------------------------ |
|---|
| 108 | ; We find common points between u and v |
|---|
| 109 | ;------------------------------------------------------------ |
|---|
| 110 | indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1] |
|---|
| 111 | indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1] |
|---|
| 112 | indicex = inter(indicexu, indicexv) |
|---|
| 113 | indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1] |
|---|
| 114 | indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1] |
|---|
| 115 | indicey = inter(indiceyu, indiceyv) |
|---|
| 116 | nx = n_elements(indicex) |
|---|
| 117 | ny = n_elements(indicey) |
|---|
| 118 | ;---------------------------------------------------------------------------- |
|---|
| 119 | vargrid = 'T' |
|---|
| 120 | varname = 'norm' |
|---|
| 121 | if n_elements(valmask) EQ 0 THEN valmask = 1e20 |
|---|
| 122 | firstxt = indicex[0] & lastxt = indicex[0]+nx-1 & nxt = nx |
|---|
| 123 | firstyt = indicey[0] & lastyt = indicey[0]+ny-1 & nyt = ny |
|---|
| 124 | ;---------------------------------------------------------------------------- |
|---|
| 125 | case 1 of |
|---|
| 126 | ;---------------------------------------------------------------------------- |
|---|
| 127 | ;---------------------------------------------------------------------------- |
|---|
| 128 | ;xyz |
|---|
| 129 | ;---------------------------------------------------------------------------- |
|---|
| 130 | ;---------------------------------------------------------------------------- |
|---|
| 131 | (size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN |
|---|
| 132 | ;---------------------------------------------------------------------------- |
|---|
| 133 | indice2d = lindgen(jpi, jpj) |
|---|
| 134 | indice2d = indice2d[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1] |
|---|
| 135 | indice3d = lindgen(jpi, jpj, jpk) |
|---|
| 136 | indice3d = indice3d[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt] |
|---|
| 137 | ;------------------------------------------------------------ |
|---|
| 138 | ; extraction of u and v on the appropriated domain |
|---|
| 139 | ;------------------------------------------------------------ |
|---|
| 140 | case 1 of |
|---|
| 141 | (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ |
|---|
| 142 | (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN |
|---|
| 143 | case (size(u))[3] OF |
|---|
| 144 | nzt:BEGIN |
|---|
| 145 | if nxu NE nx then $ |
|---|
| 146 | if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] |
|---|
| 147 | IF nxv NE nx THEN $ |
|---|
| 148 | if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] |
|---|
| 149 | IF nyu NE ny THEN $ |
|---|
| 150 | if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] |
|---|
| 151 | IF nyv NE ny THEN $ |
|---|
| 152 | if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] |
|---|
| 153 | end |
|---|
| 154 | jpk:BEGIN |
|---|
| 155 | if nxu NE nx then $ |
|---|
| 156 | if indicex[0] EQ firstxu then u = u[0:nx-1, *, firstzt:lastzt] ELSE u = u[1: nx, *, firstzt:lastzt] |
|---|
| 157 | IF nxv NE nx THEN $ |
|---|
| 158 | if indicex[0] EQ firstxv then v = v[0:nx-1, *, firstzt:lastzt] ELSE v = v[1: nx, *, firstzt:lastzt] |
|---|
| 159 | IF nyu NE ny THEN $ |
|---|
| 160 | if indicey[0] EQ firstyu then u = u[*, 0:ny-1, firstzt:lastzt] ELSE u = u[*, 1: ny, firstzt:lastzt] |
|---|
| 161 | IF nyv NE ny THEN $ |
|---|
| 162 | if indicey[0] EQ firstyv then v = v[*, 0:ny-1, firstzt:lastzt] ELSE v = v[*, 1: ny, firstzt:lastzt] |
|---|
| 163 | end |
|---|
| 164 | ELSE: return, report(['the third dimension of u (' + strtrim((size(u))[3], 1) $ |
|---|
| 165 | +') must be equal to nzt (' + strtrim(nzt, 1) + ') or jpk ('+strtrim(jpk, 1)+')']) |
|---|
| 166 | endcase |
|---|
| 167 | END |
|---|
| 168 | (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND (size(u))[3] EQ jpk AND $ |
|---|
| 169 | (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj AND (size(u))[3] EQ jpk :BEGIN |
|---|
| 170 | u = u[indice3d] |
|---|
| 171 | v = v[indice3d] |
|---|
| 172 | END |
|---|
| 173 | ELSE: return $ |
|---|
| 174 | , report(['We cannot find re-cut on common indexes (see restriction un the header.).' $ |
|---|
| 175 | , 'To avoid this problem, read the full domain' $ |
|---|
| 176 | , 'or use the keyword /memeindice in domdef when defining the zoom area']) |
|---|
| 177 | endcase |
|---|
| 178 | ;------------------------------------------------------------------ |
|---|
| 179 | ; We reshape u and v to make sure that no dimension has been erased |
|---|
| 180 | ;------------------------------------------------------------------ |
|---|
| 181 | if nzt EQ 1 then begin |
|---|
| 182 | u = reform(u, nx, ny, nzt, /over) |
|---|
| 183 | v = reform(v, nx, ny, nzt, /over) |
|---|
| 184 | endif |
|---|
| 185 | ;------------------------------------------------------------------ |
|---|
| 186 | ; construction of u and v at points T |
|---|
| 187 | ;----------------------------------------------------------- |
|---|
| 188 | a = u[0, *, *] |
|---|
| 189 | u = (u+shift(u, 1, 0, 0))/2. |
|---|
| 190 | if NOT keyword_set(key_periodic) OR nx NE jpi then u[0, *, *] = a |
|---|
| 191 | a = v[*, 0, *] |
|---|
| 192 | v = (v+shift(v, 0, 1, 0))/2. |
|---|
| 193 | if NOT keyword_set(key_periodic) OR nx NE jpi then v[*, 0, *] = a |
|---|
| 194 | ;---------------------------------------------------------------------------- |
|---|
| 195 | ; attribution of the mask and of longitude and latitude arrays |
|---|
| 196 | ;---------------------------------------------------------------------------- |
|---|
| 197 | mask = tmask[indice3d] |
|---|
| 198 | if nzt EQ 1 then mask = reform(mask, nx, ny, nzt, /over) |
|---|
| 199 | ;----------------------------------------------------------- |
|---|
| 200 | if n_elements(valmask) EQ 0 THEN valmask = 1e20 |
|---|
| 201 | landu = where(u GE valmask/10) |
|---|
| 202 | if landu[0] NE -1 then u[landu] = 0 |
|---|
| 203 | landv = where(v GE valmask/10) |
|---|
| 204 | if landv[0] NE -1 then v[landv] = 0 |
|---|
| 205 | res = sqrt(u^2+v^2) |
|---|
| 206 | if NOT keyword_set(key_periodic) OR nx NE jpi then res[0, *, *] = !values.f_nan |
|---|
| 207 | res[*, 0, *] = !values.f_nan |
|---|
| 208 | mask = where(mask eq 0) |
|---|
| 209 | IF mask[0] NE -1 THEN res[mask] = valmask |
|---|
| 210 | ; All kind of average |
|---|
| 211 | domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0], (gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme |
|---|
| 212 | if keyword_set(direc) then res = moyenne(res, direc, /nan, boxzoom = boxzoom, /nodomdef) |
|---|
| 213 | ;----------------------------------------------------------- |
|---|
| 214 | END |
|---|
| 215 | ;---------------------------------------------------------------------------- |
|---|
| 216 | ;---------------------------------------------------------------------------- |
|---|
| 217 | ;xyt |
|---|
| 218 | ;---------------------------------------------------------------------------- |
|---|
| 219 | ;---------------------------------------------------------------------------- |
|---|
| 220 | date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN |
|---|
| 221 | indice2d = lindgen(jpi, jpj) |
|---|
| 222 | indice2d = indice2d[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1] |
|---|
| 223 | ;------------------------------------------------------------ |
|---|
| 224 | ; extraction of u and v on the appropriated domain |
|---|
| 225 | ;------------------------------------------------------------ |
|---|
| 226 | case 1 of |
|---|
| 227 | (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ |
|---|
| 228 | (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN |
|---|
| 229 | if nxu NE nx then $ |
|---|
| 230 | if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] |
|---|
| 231 | IF nxv NE nx THEN $ |
|---|
| 232 | if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] |
|---|
| 233 | IF nyu NE ny THEN $ |
|---|
| 234 | if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] |
|---|
| 235 | IF nyv NE ny THEN $ |
|---|
| 236 | if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] |
|---|
| 237 | END |
|---|
| 238 | (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ |
|---|
| 239 | (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN |
|---|
| 240 | u = u[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *] |
|---|
| 241 | v = v[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *] |
|---|
| 242 | END |
|---|
| 243 | ELSE:return $ |
|---|
| 244 | , report(['We cannot find re-cut on common indexes (see restriction un the header.).' $ |
|---|
| 245 | , 'To avoid this problem, read the full domain' $ |
|---|
| 246 | , 'or use the keyword /memeindice in domdef when defining the zoom area']) |
|---|
| 247 | endcase |
|---|
| 248 | ;------------------------------------------------------------------ |
|---|
| 249 | ; construction of u and v at points T |
|---|
| 250 | ;----------------------------------------------------------- |
|---|
| 251 | a = u[0, *, *] |
|---|
| 252 | u = (u+shift(u, 1, 0, 0))/2. |
|---|
| 253 | if NOT keyword_set(key_periodic) OR nx NE jpi then u[0, *, *] = a |
|---|
| 254 | a = v[*, 0, *] |
|---|
| 255 | v = (v+shift(v, 0, 1, 0))/2. |
|---|
| 256 | if NOT keyword_set(key_periodic) OR nx NE jpi then v[*, 0, *] = a |
|---|
| 257 | ;---------------------------------------------------------------------------- |
|---|
| 258 | ; attribution of the mask and of longitude and latitude arrays. |
|---|
| 259 | ; We recover the complete grid to establish a big mask extent in the four |
|---|
| 260 | ; direction to cover pointsfor which a land point has been |
|---|
| 261 | ; considerated (make a small drawing) |
|---|
| 262 | ;---------------------------------------------------------------------------- |
|---|
| 263 | mask = tmask[indice2d+jpi*jpj*firstzt] |
|---|
| 264 | if ny EQ 1 then mask = reform(mask, nx, ny, /over) |
|---|
| 265 | ;----------------------------------------------------------- |
|---|
| 266 | ; construction of land containing all points to mask |
|---|
| 267 | ;----------------------------------------------------------- |
|---|
| 268 | if n_elements(valmask) EQ 0 THEN valmask = 1e20 |
|---|
| 269 | landu = where(u GE valmask/10) |
|---|
| 270 | if landu[0] NE -1 then u[landu] = 0 |
|---|
| 271 | landv = where(v GE valmask/10) |
|---|
| 272 | if landv[0] NE -1 then v[landv] = 0 |
|---|
| 273 | res = sqrt(u^2+v^2) |
|---|
| 274 | if NOT keyword_set(key_periodic) OR nx NE jpi then res[0, *, *] = !values.f_nan |
|---|
| 275 | res[*, 0, *] = !values.f_nan |
|---|
| 276 | mask = where(mask eq 0) |
|---|
| 277 | IF mask[0] NE -1 THEN BEGIN |
|---|
| 278 | coeftps = lindgen(jpt)*nx*ny |
|---|
| 279 | coeftps = replicate(1, n_elements(mask))#coeftps |
|---|
| 280 | mask = (temporary(mask))[*]#replicate(1, jpt) |
|---|
| 281 | mask = temporary(mask[*]) + temporary(coeftps[*]) |
|---|
| 282 | res[temporary(mask)] = valmask |
|---|
| 283 | ENDIF |
|---|
| 284 | ; moyennes en tous genres |
|---|
| 285 | domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0], (gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme |
|---|
| 286 | if keyword_set(direc) then res = grossemoyenne(res, direc, /nan, boxzoom = boxzoom, /nodomdef) |
|---|
| 287 | END |
|---|
| 288 | ;---------------------------------------------------------------------------- |
|---|
| 289 | ;---------------------------------------------------------------------------- |
|---|
| 290 | ;xyzt |
|---|
| 291 | ;---------------------------------------------------------------------------- |
|---|
| 292 | ;---------------------------------------------------------------------------- |
|---|
| 293 | date1 NE date2 AND (size(u))[0] EQ 4:BEGIN |
|---|
| 294 | indice2d = lindgen(jpi, jpj) |
|---|
| 295 | indice2d = indice2d[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1] |
|---|
| 296 | indice3d = lindgen(jpi, jpj, jpk) |
|---|
| 297 | indice3d = indice3d[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt] |
|---|
| 298 | ;------------------------------------------------------------ |
|---|
| 299 | ; extraction of u and v on the appropriated domain |
|---|
| 300 | ;------------------------------------------------------------ |
|---|
| 301 | case 1 of |
|---|
| 302 | (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ |
|---|
| 303 | (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN |
|---|
| 304 | case (size(u))[3] OF |
|---|
| 305 | nzt:BEGIN |
|---|
| 306 | if nxu NE nx then $ |
|---|
| 307 | if indicex[0] EQ firstxu then u = u[0:nx-1, *, *, *] ELSE u = u[1: nx, *, *, *] |
|---|
| 308 | IF nxv NE nx THEN $ |
|---|
| 309 | if indicex[0] EQ firstxv then v = v[0:nx-1, *, *, *] ELSE v = v[1: nx, *, *, *] |
|---|
| 310 | IF nyu NE ny THEN $ |
|---|
| 311 | if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *, *] ELSE u = u[*, 1: ny, *, *] |
|---|
| 312 | IF nyv NE ny THEN $ |
|---|
| 313 | if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *, *] ELSE v = v[*, 1: ny, *, *] |
|---|
| 314 | end |
|---|
| 315 | jpk:BEGIN |
|---|
| 316 | if nxu NE nx then $ |
|---|
| 317 | if indicex[0] EQ firstxu then u = u[0:nx-1, *, firstzt:lastzt, *] ELSE u = u[1: nx, *, firstzt:lastzt, *] |
|---|
| 318 | IF nxv NE nx THEN $ |
|---|
| 319 | if indicex[0] EQ firstxv then v = v[0:nx-1, *, firstzt:lastzt, *] ELSE v = v[1: nx, *, firstzt:lastzt, *] |
|---|
| 320 | IF nyu NE ny THEN $ |
|---|
| 321 | if indicey[0] EQ firstyu then u = u[*, 0:ny-1, firstzt:lastzt, *] ELSE u = u[*, 1: ny, firstzt:lastzt, *] |
|---|
| 322 | IF nyv NE ny THEN $ |
|---|
| 323 | if indicey[0] EQ firstyv then v = v[*, 0:ny-1, firstzt:lastzt, *] ELSE v = v[*, 1: ny, firstzt:lastzt, *] |
|---|
| 324 | end |
|---|
| 325 | ELSE: report, (['the third dimension of u (' + strtrim((size(u))[3], 1) $ |
|---|
| 326 | +') must be equal to nzt (' + strtrim(nzt, 1) + ') or jpk ('+strtrim(jpk, 1)+')']) |
|---|
| 327 | endcase |
|---|
| 328 | END |
|---|
| 329 | (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND (size(u))[3] EQ jpk AND $ |
|---|
| 330 | (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj AND (size(u))[3] EQ jpk :BEGIN |
|---|
| 331 | u = u[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt, *] |
|---|
| 332 | v = v[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt, *] |
|---|
| 333 | END |
|---|
| 334 | ELSE: return $ |
|---|
| 335 | , report(['We cannot find re-cut on common indexes (see restriction un the header.).' $ |
|---|
| 336 | , 'To avoid this problem, read the full domain' $ |
|---|
| 337 | , 'or use the keyword /memeindice in domdef when defining the zoom area']) |
|---|
| 338 | endcase |
|---|
| 339 | ;------------------------------------------------------------------ |
|---|
| 340 | ; construction of u and v at points T |
|---|
| 341 | ;----------------------------------------------------------- |
|---|
| 342 | a = u[0, *, *, *] |
|---|
| 343 | u = (u+shift(u, 1, 0, 0, 0))/2. |
|---|
| 344 | if NOT keyword_set(key_periodic) OR nx NE jpi then u[0, *, *, *] = a |
|---|
| 345 | a = v[*, 0, *, *] |
|---|
| 346 | v = (v+shift(v, 0, 1, 0, 0))/2. |
|---|
| 347 | if NOT keyword_set(key_periodic) OR nx NE jpi then v[*, 0, *, *] = a |
|---|
| 348 | ;---------------------------------------------------------------------------- |
|---|
| 349 | ; attribution of the mask and of longitude and latitude arrays |
|---|
| 350 | ;---------------------------------------------------------------------------- |
|---|
| 351 | mask = tmask[indice3d] |
|---|
| 352 | if nzt EQ 1 then mask = reform(mask, nx, ny, nzt, /over) |
|---|
| 353 | ;----------------------------------------------------------- |
|---|
| 354 | if n_elements(valmask) EQ 0 THEN valmask = 1e20 |
|---|
| 355 | landu = where(u GE valmask/10) |
|---|
| 356 | if landu[0] NE -1 then u[landu] = 0 |
|---|
| 357 | landv = where(v GE valmask/10) |
|---|
| 358 | if landv[0] NE -1 then v[landv] = 0 |
|---|
| 359 | res = sqrt(u^2+v^2) |
|---|
| 360 | if NOT keyword_set(key_periodic) OR nx NE jpi then res[0, *, *, *] = !values.f_nan |
|---|
| 361 | res[*, 0, *, *] = !values.f_nan |
|---|
| 362 | mask = where(mask eq 0) |
|---|
| 363 | IF mask[0] NE -1 THEN BEGIN |
|---|
| 364 | coeftps = lindgen(jpt)*nx*ny*nzt |
|---|
| 365 | coeftps = replicate(1, n_elements(mask))#coeftps |
|---|
| 366 | mask = (temporary(mask))[*]#replicate(1, jpt) |
|---|
| 367 | mask = temporary(mask[*]) + temporary(coeftps[*]) |
|---|
| 368 | res[temporary(mask)] = valmask |
|---|
| 369 | ENDIF |
|---|
| 370 | ; All kind of average |
|---|
| 371 | domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0], (gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme |
|---|
| 372 | if keyword_set(direc) then res = grossemoyenne(res, direc, /nan, boxzoom = boxzoom, /nodomdef) |
|---|
| 373 | END |
|---|
| 374 | ;---------------------------------------------------------------------------- |
|---|
| 375 | ;---------------------------------------------------------------------------- |
|---|
| 376 | ;xy |
|---|
| 377 | ;---------------------------------------------------------------------------- |
|---|
| 378 | ;---------------------------------------------------------------------------- |
|---|
| 379 | ELSE:BEGIN ;xy |
|---|
| 380 | indice2d = lindgen(jpi, jpj) |
|---|
| 381 | indice2d = indice2d[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1] |
|---|
| 382 | ;------------------------------------------------------------ |
|---|
| 383 | ; extraction of u and v on the appropriated domain |
|---|
| 384 | ;------------------------------------------------------------ |
|---|
| 385 | case 1 of |
|---|
| 386 | (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ |
|---|
| 387 | (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN |
|---|
| 388 | if nxu NE nx then $ |
|---|
| 389 | if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] |
|---|
| 390 | IF nxv NE nx THEN $ |
|---|
| 391 | if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] |
|---|
| 392 | IF nyu NE ny THEN $ |
|---|
| 393 | if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] |
|---|
| 394 | IF nyv NE ny THEN $ |
|---|
| 395 | if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] |
|---|
| 396 | END |
|---|
| 397 | (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ |
|---|
| 398 | (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN |
|---|
| 399 | u = u[indice2d] |
|---|
| 400 | v = v[indice2d] |
|---|
| 401 | END |
|---|
| 402 | ELSE:return $ |
|---|
| 403 | , report(['We cannot find re-cut on common indexes (see restriction un the header.).' $ |
|---|
| 404 | , 'To avoid this problem, read the full domain' $ |
|---|
| 405 | , 'or use the keyword /memeindice in domdef when defining the zoom area']) |
|---|
| 406 | endcase |
|---|
| 407 | ;------------------------------------------------------------------ |
|---|
| 408 | ; We reshape u and v to make sure that no dimension has been erased |
|---|
| 409 | ;------------------------------------------------------------------ |
|---|
| 410 | if ny EQ 1 then begin |
|---|
| 411 | u = reform(u, nx, ny, /over) |
|---|
| 412 | v = reform(v, nx, ny, /over) |
|---|
| 413 | endif |
|---|
| 414 | ;------------------------------------------------------------------ |
|---|
| 415 | ; construction of u and v at points T |
|---|
| 416 | ;----------------------------------------------------------- |
|---|
| 417 | a = u[0, *] |
|---|
| 418 | u = (u+shift(u, 1, 0))/2. |
|---|
| 419 | if NOT keyword_set(key_periodic) OR nx NE jpi then u[0, *] = a |
|---|
| 420 | a = v[*, 0] |
|---|
| 421 | v = (v+shift(v, 0, 1))/2. |
|---|
| 422 | if NOT keyword_set(key_periodic) OR nx NE jpi then v[*, 0] = a |
|---|
| 423 | ;---------------------------------------------------------------------------- |
|---|
| 424 | ; attribution of the mask and of longitude and latitude arrays. |
|---|
| 425 | ; We recover the complete grid to establish a big mask extent in the four |
|---|
| 426 | ; direction to cover pointsfor which a land point has been |
|---|
| 427 | ; considerated (make a small drawing) |
|---|
| 428 | ;---------------------------------------------------------------------------- |
|---|
| 429 | mask = tmask[indice2d+jpi*jpj*firstzt] |
|---|
| 430 | if nyt EQ 1 THEN mask = reform(mask, nx, ny, /over) |
|---|
| 431 | ;----------------------------------------------------------- |
|---|
| 432 | ; construction of land containing all points to mask |
|---|
| 433 | ;----------------------------------------------------------- |
|---|
| 434 | if n_elements(valmask) EQ 0 THEN valmask = 1e20 |
|---|
| 435 | landu = where(u GE valmask/10) |
|---|
| 436 | if landu[0] NE -1 then u[landu] = 0 |
|---|
| 437 | landv = where(v GE valmask/10) |
|---|
| 438 | if landv[0] NE -1 then v[landv] = 0 |
|---|
| 439 | res = sqrt(u^2+v^2) |
|---|
| 440 | if NOT keyword_set(key_periodic) OR nx NE jpi then res[0, *] = !values.f_nan |
|---|
| 441 | res[*, 0] = !values.f_nan |
|---|
| 442 | mask = where(mask eq 0) |
|---|
| 443 | IF mask[0] NE -1 THEN res[mask] = valmask |
|---|
| 444 | ; All kind of average |
|---|
| 445 | domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0], (gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme |
|---|
| 446 | if keyword_set(direc) then res = moyenne(res, direc, /nan, boxzoom = boxzoom, /nodomdef) |
|---|
| 447 | END |
|---|
| 448 | ;---------------------------------------------------------------------------- |
|---|
| 449 | endcase |
|---|
| 450 | ;------------------------------------------------------------ |
|---|
| 451 | if keyword_set(key_performance) THEN print, 'time norm', systime(1)-time1 |
|---|
| 452 | return, res |
|---|
| 453 | end |
|---|
