Ignore:
Timestamp:
12/03/07 15:12:28 (17 years ago)
Author:
smasson
Message:

bugfix + cleaning of Computation routines

File:
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/SRC/Computation/norm.pro

    r312 r314  
    22; 
    33; @file_comments 
    4 ; calculate the norm of a field of vectors, then make a possible average. 
    5 ;   Comment 1: The field of vector can be, 2d:xy, 3d: xyz or xyt, 
    6 ; 4d: xyzt 
    7 ;   Comment 2: 
    8 ; The calculation of the norm is made before the possible spatial or 
    9 ; temporal average because the average of the norm is not equal to the 
    10 ; norm of averages 
     4; calculate the norm of vectors field located on Arakawa C-grid 
    115; 
    126; @categories 
    137; Calculation 
    148; 
    15 ; @param COMPOSANTEU {in}{required} 
    16 ; an 2d, 3d or 4d array 
    17 ; 
    18 ; @param COMPOSANTEV {in}{required} 
    19 ; an 2d, 3d or 4d array 
    20 ; 
    21 ; @keyword BOXZOOM 
    22 ; boxzoom on which do the average (by default the domain selected 
    23 ; by the last <pro>domdef</pro> done) 
     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. 
    2420; 
    2521; @keyword DIREC 
    … …  
    2723;       'xzt' 'yzt' 'xyzt' Direction on which do averages 
    2824; 
     25; @keyword _EXTRA 
     26; Used to declare that this routine accepts inherited keywords 
     27; 
    2928; @returns 
    30 ; Array to trace with <pro>plt</pro>, <pro>pltz</pro> or <pro>pltt</pro>. 
     29; A 2D (xy), 3D (xyz or yt), 4D (xyzt) Array 
    3130; 
    3231; @uses 
    33 ; common.pro 
     32; @cm_4mesh 
     33; @cm_4data 
     34; @cm_4cal 
    3435; 
    3536; @restrictions 
    3637; The norm is calculated on points T. To do this calculation, we average 
    37 ; field U and V on points T before calculate the norme. At the edge of 
     38; field U and V on points T before calculate the norm. At the edge of 
    3839; coast and of domain, we can not calculate fields U and V at points T, 
    3940; that is why these points are at value !values.f_nan. 
    … …  
    4647; To know what type of array we work with, we  test its size and dates 
    4748; gave by time[0] and time[jpt-1] to know if thee is a temporal dimension. 
    48 ; Before to start norme, make sure that time and jpt are defined how 
     49; Before to start norm, make sure that time and jpt are defined how 
    4950; they have to! 
    5051; 
    5152; @examples 
    52 ; To calculate the average of the norme of streams on all the domain 
     53; To calculate the average of the norm of streams on all the domain 
    5354; between 0 and 50: 
    54 ;      IDL> res=norme(un,vn,boxzoom=[0,50],dir='xyz') 
     55;      IDL> domdef, 0, 50 
     56;      IDL> res = norm(un, vn, dir = 'xyz') 
    5557; 
    5658; @history 
    5759; Sebastien Masson (smasson\@lodyc.jussieu.fr) 
    5860;                       9/6/1999 
    59 ; 
    6061; @version 
    6162; $Id$ 
    … …  
    6364;- 
    6465; 
    65 FUNCTION norme, composanteu, composantev, BOXZOOM = boxzoom, DIREC = direc, _EXTRA = ex 
     66FUNCTION norm, uu, vv, DIREC = direc, _EXTRA = ex 
    6667; 
    6768  compile_opt idl2, strictarrsubs 
    … …  
    7576  ENDIF 
    7677;--------------------------------------------------------- 
    77    tempsun = systime(1)         ; To key_performance 
    78 ; 
    79    IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 
    80      return, report(['This version of norme is based on Arakawa C-grid.' $ 
     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.' $ 
    8182                     , 'U and V grids must therefore be defined']) 
    82 ; 
    83 ;------------------------------------------------------------ 
    84   if keyword_set(boxzoom) then BEGIN 
    85     Case 1 Of 
    86       N_Elements(Boxzoom) Eq 1:bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]] 
    87       N_Elements(Boxzoom) Eq 2:bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]] 
    88       N_Elements(Boxzoom) Eq 4:bte = [Boxzoom, vert1, vert2] 
    89       N_Elements(Boxzoom) Eq 5:bte = [Boxzoom[0:3], 0, Boxzoom[4]] 
    90       N_Elements(Boxzoom) Eq 6:bte = Boxzoom 
    91       Else: return, report('Mauvaise Definition de Boxzoom') 
    92     ENDCASE 
    93     domdef, boxzoom 
    94   ENDIF 
    95 ;------------------------------------------------------------ 
    96    if NOT keyword_set(direc) then direc = 0 
     83;------------------------------------------------------------ 
     84  if NOT keyword_set(direc) then direc = 0 
    9785; construction of u and v at points T 
    98    u = litchamp(composanteu) 
    99    v = litchamp(composantev) 
    100    date1 = time[0] 
    101    if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1] 
     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] 
    10290 
    103    if (size(u))[0] NE (size(v))[0] then return,  -1 
    104  
    105    vargrid='T' 
    106    varname = 'norme ' 
    107    valmask = 1e20 
    108 ; 
    109    grilleu = litchamp(composanteu, /grid) 
    110    if grilleu EQ '' then grilleu = 'U' 
    111    grillev = litchamp(composantev, /grid) 
    112    if grillev EQ '' then grillev = 'V' 
    113    IF grilleu EQ 'V' AND grillev EQ 'U' THEN inverse = 1 
    114    IF grilleu EQ 'T' AND grillev EQ 'T' THEN BEGIN 
    115       interpolle  = 0 
    116       return, report('cas non code mais facile a faire!') 
    117    ENDIF ELSE interpolle = 1 
    118    if keyword_set(inverse) then begin 
    119       rien = u 
    120       u = v 
    121       v = rien 
    122    endif 
    123  
    124  
     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 
    125107;------------------------------------------------------------ 
    126108; We find common points between u and v 
    127109;------------------------------------------------------------ 
    128    indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1] 
    129    indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1] 
    130    indicex = inter(indicexu, indicexv) 
    131    indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1] 
    132    indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1] 
    133    indicey = inter(indiceyu, indiceyv) 
    134    nx = n_elements(indicex) 
    135    ny = n_elements(indicey) 
    136 ;---------------------------------------------------------------------------- 
    137    case 1 of 
     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 
    138126;---------------------------------------------------------------------------- 
    139127;---------------------------------------------------------------------------- 
    … …  
    141129;---------------------------------------------------------------------------- 
    142130;---------------------------------------------------------------------------- 
    143       (size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN 
    144 ;---------------------------------------------------------------------------- 
    145          indice2d = lindgen(jpi, jpj) 
    146          indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1] 
    147          indice3d = lindgen(jpi, jpj, jpk) 
    148          indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt] 
     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] 
    149137;------------------------------------------------------------ 
    150138; extraction of u and v on the appropriated domain 
    151139;------------------------------------------------------------ 
    152          case 1 of 
    153             (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 
    154              (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
    155                case (size(u))[3] OF 
    156                   nzt:BEGIN 
    157                      if nxu NE nx then $ 
    158                       if indicex[0] EQ firstxu then u = u[0:nx-1,*,*] ELSE u = u[1: nx, *,*] 
    159                      IF nxv NE nx THEN $ 
    160                       if indicex[0] EQ firstxv then v = v[0:nx-1,*,*] ELSE v = v[1: nx, *,*] 
    161                      IF nyu NE ny THEN $ 
    162                       if indicey[0] EQ firstyu then u = u[*,0:ny-1,*] ELSE u = u[*, 1: ny,*] 
    163                      IF nyv NE ny THEN $ 
    164                       if indicey[0] EQ firstyv then v = v[*,0:ny-1,*] ELSE v = v[*, 1: ny,*] 
    165                   end 
    166                   jpk:BEGIN 
    167                      if nxu NE nx then $ 
    168                       if indicex[0] EQ firstxu then u = u[0:nx-1, *,firstzt:lastzt] ELSE u = u[1: nx, *,firstzt:lastzt] 
    169                      IF nxv NE nx THEN $ 
    170                       if indicex[0] EQ firstxv then v = v[0:nx-1, *,firstzt:lastzt] ELSE v = v[1: nx, *,firstzt:lastzt] 
    171                      IF nyu NE ny THEN $ 
    172                       if indicey[0] EQ firstyu then u = u[*, 0:ny-1,firstzt:lastzt] ELSE u = u[*, 1: ny,firstzt:lastzt] 
    173                      IF nyv NE ny THEN $ 
    174                       if indicey[0] EQ firstyv then v = v[*, 0:ny-1,firstzt:lastzt] ELSE v = v[*, 1: ny,firstzt:lastzt] 
    175                   end 
    176                   ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs') 
    177                endcase 
    178             END 
    179             (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND (size(u))[3] EQ jpk AND $ 
    180              (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj AND (size(u))[3] EQ jpk :BEGIN 
    181                u = u[indice3d] 
    182                v = v[indice3d] 
    183             END 
    184             ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs') 
    185          endcase 
     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 
    186178;------------------------------------------------------------------ 
    187179; We reshape u and v to make sure that no dimension has been erased 
    188180;------------------------------------------------------------------ 
    189          if nzt EQ 1 then begin 
    190             u = reform(u, nx, ny, nzt, /over) 
    191             v = reform(v, nx, ny, nzt, /over) 
    192          endif 
     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 
    193185;------------------------------------------------------------------ 
    194186; construction of u and v at points T 
    195187;----------------------------------------------------------- 
    196          a=u[0,*,*] 
    197          u=(u+shift(u,1,0,0))/2. 
    198          if NOT keyword_set(key_periodic) OR nx NE jpi then u[0,*,*]=a 
    199          a=v[*,0,*] 
    200          v=(v+shift(v,0,1,0))/2. 
    201          if NOT keyword_set(key_periodic) OR nx NE jpi then v[*,0,*]=a 
     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 
    202194;---------------------------------------------------------------------------- 
    203195; attribution of the mask and of longitude and latitude arrays 
    204196;---------------------------------------------------------------------------- 
    205          mask = tmask[indice3d] 
    206          if nzt EQ 1 then mask = reform(mask, nx, ny, nzt, /over) 
    207 ;----------------------------------------------------------- 
    208          if n_elements(valmask) EQ 0 THEN valmask = 1e20 
    209          landu = where(u GE valmask/10) 
    210          if landu[0] NE -1 then u[landu] = 0 
    211          landv = where(v GE valmask/10) 
    212          if landv[0] NE -1 then v[landv] = 0 
    213          res=sqrt(u^2+v^2) 
    214          if NOT keyword_set(key_periodic) OR nx NE jpi then res[0,*, *]=!values.f_nan 
    215          res[*,0, *]=!values.f_nan 
    216          mask = where(mask eq 0) 
    217          IF mask[0] NE -1 THEN res[mask] = valmask 
     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 
    218210; All kind of average 
    219          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme 
    220          if keyword_set(direc) then res = moyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef) 
    221 ;----------------------------------------------------------- 
    222       END 
     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 
    223215;---------------------------------------------------------------------------- 
    224216;---------------------------------------------------------------------------- 
    … …  
    226218;---------------------------------------------------------------------------- 
    227219;---------------------------------------------------------------------------- 
    228       date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN 
    229          indice2d = lindgen(jpi, jpj) 
    230          indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1] 
     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] 
    231223;------------------------------------------------------------ 
    232224; extraction of u and v on the appropriated domain 
    233225;------------------------------------------------------------ 
    234          case 1 of 
    235             (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 
    236              (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
    237                if nxu NE nx then $ 
    238                 if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
    239                IF nxv NE nx THEN $ 
    240                 if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
    241                IF nyu NE ny THEN $ 
    242                 if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
    243                IF nyv NE ny THEN $ 
    244                 if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
    245             END 
    246             (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 
    247              (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
    248                u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 
    249                v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 
    250             END 
    251             ELSE:return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs') 
    252          endcase 
     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 
    253248;------------------------------------------------------------------ 
    254249; construction of u and v at points T 
    255250;----------------------------------------------------------- 
    256          a=u[0,*,*] 
    257          u=(u+shift(u,1,0,0))/2. 
    258          if NOT keyword_set(key_periodic) OR nx NE jpi then u[0,*,*]=a 
    259          a=v[*,0,*] 
    260          v=(v+shift(v,0,1,0))/2. 
    261          if NOT keyword_set(key_periodic) OR nx NE jpi then v[*,0,*]=a 
     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 
    262257;---------------------------------------------------------------------------- 
    263258; attribution of the mask and of longitude and latitude arrays. 
    … …  
    266261; considerated (make a small drawing) 
    267262;---------------------------------------------------------------------------- 
    268          mask = tmask[indice2d+jpi*jpj*firstzt] 
    269          if ny EQ 1 then mask = reform(mask, nx, ny, /over) 
     263      mask = tmask[indice2d+jpi*jpj*firstzt] 
     264      if ny EQ 1 then mask = reform(mask, nx, ny, /over) 
    270265;----------------------------------------------------------- 
    271266; construction of land containing all points to mask 
    272267;----------------------------------------------------------- 
    273          if n_elements(valmask) EQ 0 THEN valmask = 1e20 
    274          landu = where(u GE valmask/10) 
    275          if landu[0] NE -1 then u[landu] = 0 
    276          landv = where(v GE valmask/10) 
    277          if landv[0] NE -1 then v[landv] = 0 
    278          res=sqrt(u^2+v^2) 
    279          if NOT keyword_set(key_periodic) OR nx NE jpi then res[0,*, *]=!values.f_nan 
    280          res[*,0, *]=!values.f_nan 
    281          mask = where(mask eq 0) 
    282          IF mask[0] NE -1 THEN BEGIN 
    283             coeftps = lindgen(jpt)*nx*ny 
    284             coeftps = replicate(1, n_elements(mask))#coeftps 
    285             mask = (temporary(mask))[*]#replicate(1, jpt) 
    286             mask =temporary(mask[*]) + temporary(coeftps[*]) 
    287             res[temporary(mask)] = valmask 
    288          ENDIF 
     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 
    289284; moyennes en tous genres 
    290          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme 
    291          if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef) 
    292       END 
     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 
    293288;---------------------------------------------------------------------------- 
    294289;---------------------------------------------------------------------------- 
    … …  
    296291;---------------------------------------------------------------------------- 
    297292;---------------------------------------------------------------------------- 
    298       date1 NE date2 AND (size(u))[0] EQ 4:BEGIN 
    299          indice2d = lindgen(jpi, jpj) 
    300          indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1] 
    301          indice3d = lindgen(jpi, jpj, jpk) 
    302          indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt] 
     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] 
    303298;------------------------------------------------------------ 
    304299; extraction of u and v on the appropriated domain 
    305300;------------------------------------------------------------ 
    306          case 1 of 
    307             (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 
    308              (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
    309                case (size(u))[3] OF 
    310                   nzt:BEGIN 
    311                      if nxu NE nx then $ 
    312                       if indicex[0] EQ firstxu then u = u[0:nx-1,*,*,*] ELSE u = u[1: nx, *,*,*] 
    313                      IF nxv NE nx THEN $ 
    314                       if indicex[0] EQ firstxv then v = v[0:nx-1,*,*,*] ELSE v = v[1: nx, *,*,*] 
    315                      IF nyu NE ny THEN $ 
    316                       if indicey[0] EQ firstyu then u = u[*,0:ny-1,*,*] ELSE u = u[*, 1: ny,*,*] 
    317                      IF nyv NE ny THEN $ 
    318                       if indicey[0] EQ firstyv then v = v[*,0:ny-1,*,*] ELSE v = v[*, 1: ny,*,*] 
    319                   end 
    320                   jpk:BEGIN 
    321                      if nxu NE nx then $ 
    322                       if indicex[0] EQ firstxu then u = u[0:nx-1, *,firstzt:lastzt,*] ELSE u = u[1: nx, *,firstzt:lastzt,*] 
    323                      IF nxv NE nx THEN $ 
    324                       if indicex[0] EQ firstxv then v = v[0:nx-1, *,firstzt:lastzt,*] ELSE v = v[1: nx, *,firstzt:lastzt,*] 
    325                      IF nyu NE ny THEN $ 
    326                       if indicey[0] EQ firstyu then u = u[*, 0:ny-1,firstzt:lastzt,*] ELSE u = u[*, 1: ny,firstzt:lastzt,*] 
    327                      IF nyv NE ny THEN $ 
    328                       if indicey[0] EQ firstyv then v = v[*, 0:ny-1,firstzt:lastzt,*] ELSE v = v[*, 1: ny,firstzt:lastzt,*] 
    329                   end 
    330                   ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs') 
    331                endcase 
    332             END 
    333             (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND (size(u))[3] EQ jpk AND $ 
    334              (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj AND (size(u))[3] EQ jpk :BEGIN 
    335                u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt, *] 
    336                v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt, *] 
    337             END 
    338             ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs') 
    339          endcase 
     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 
    340339;------------------------------------------------------------------ 
    341340; construction of u and v at points T 
    342341;----------------------------------------------------------- 
    343          a=u[0,*,*,*] 
    344          u=(u+shift(u,1,0,0,0))/2. 
    345          if NOT keyword_set(key_periodic) OR nx NE jpi then u[0,*,*,*]=a 
    346          a=v[*,0,*,*] 
    347          v=(v+shift(v,0,1,0,0))/2. 
    348          if NOT keyword_set(key_periodic) OR nx NE jpi then v[*,0,*,*]=a 
     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 
    349348;---------------------------------------------------------------------------- 
    350349; attribution of the mask and of longitude and latitude arrays 
    351350;---------------------------------------------------------------------------- 
    352          mask = tmask[indice3d] 
    353          if nzt EQ 1 then mask = reform(mask, nx, ny, nzt, /over) 
    354 ;----------------------------------------------------------- 
    355          if n_elements(valmask) EQ 0 THEN valmask = 1e20 
    356          landu = where(u GE valmask/10) 
    357          if landu[0] NE -1 then u[landu] = 0 
    358          landv = where(v GE valmask/10) 
    359          if landv[0] NE -1 then v[landv] = 0 
    360          res=sqrt(u^2+v^2) 
    361          if NOT keyword_set(key_periodic) OR nx NE jpi then res[0,*, *, *]=!values.f_nan 
    362          res[*,0, *, *]=!values.f_nan 
    363          mask = where(mask eq 0) 
    364          IF mask[0] NE -1 THEN BEGIN 
    365             coeftps = lindgen(jpt)*nx*ny*nzt 
    366             coeftps = replicate(1, n_elements(mask))#coeftps 
    367             mask = (temporary(mask))[*]#replicate(1, jpt) 
    368             mask =temporary(mask[*]) + temporary(coeftps[*]) 
    369             res[temporary(mask)] = valmask 
    370          ENDIF 
     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 
    371370; All kind of average 
    372          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme 
    373          if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef) 
    374       END 
     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 
    375374;---------------------------------------------------------------------------- 
    376375;---------------------------------------------------------------------------- 
    … …  
    378377;---------------------------------------------------------------------------- 
    379378;---------------------------------------------------------------------------- 
    380       ELSE:BEGIN                ;xy 
    381          indice2d = lindgen(jpi, jpj) 
    382          indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1] 
     379    ELSE:BEGIN                  ;xy 
     380      indice2d = lindgen(jpi, jpj) 
     381      indice2d = indice2d[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1] 
    383382;------------------------------------------------------------ 
    384383; extraction of u and v on the appropriated domain 
    385384;------------------------------------------------------------ 
    386          case 1 of 
    387             (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 
    388              (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
    389                if nxu NE nx then $ 
    390                 if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 
    391                IF nxv NE nx THEN $ 
    392                 if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 
    393                IF nyu NE ny THEN $ 
    394                 if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 
    395                IF nyv NE ny THEN $ 
    396                 if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 
    397             END 
    398             (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 
    399              (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
    400                u = u[indice2d] 
    401                v = v[indice2d] 
    402             END 
    403             ELSE:return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs') 
    404          endcase 
     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 
    405407;------------------------------------------------------------------ 
    406408; We reshape u and v to make sure that no dimension has been erased 
    407409;------------------------------------------------------------------ 
    408          if ny EQ 1 then begin 
    409             u = reform(u, nx, ny, /over) 
    410             v = reform(v, nx, ny, /over) 
    411          endif 
     410      if ny EQ 1 then begin 
     411        u = reform(u, nx, ny, /over) 
     412        v = reform(v, nx, ny, /over) 
     413      endif 
    412414;------------------------------------------------------------------ 
    413415; construction of u and v at points T 
    414416;----------------------------------------------------------- 
    415          a=u[0,*] 
    416          u=(u+shift(u,1,0))/2. 
    417          if NOT keyword_set(key_periodic) OR nx NE jpi then u[0,*]=a 
    418          a=v[*,0] 
    419          v=(v+shift(v,0,1))/2. 
    420          if NOT keyword_set(key_periodic) OR nx NE jpi then v[*,0]=a 
     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 
    421423;---------------------------------------------------------------------------- 
    422424; attribution of the mask and of longitude and latitude arrays. 
    … …  
    425427; considerated (make a small drawing) 
    426428;---------------------------------------------------------------------------- 
    427          mask = tmask[indice2d+jpi*jpj*firstzt] 
    428          if nyt EQ 1 THEN mask = reform(mask, nx, ny, /over) 
     429      mask = tmask[indice2d+jpi*jpj*firstzt] 
     430      if nyt EQ 1 THEN mask = reform(mask, nx, ny, /over) 
    429431;----------------------------------------------------------- 
    430432; construction of land containing all points to mask 
    431433;----------------------------------------------------------- 
    432          if n_elements(valmask) EQ 0 THEN valmask = 1e20 
    433          landu = where(u GE valmask/10) 
    434          if landu[0] NE -1 then u[landu] = 0 
    435          landv = where(v GE valmask/10) 
    436          if landv[0] NE -1 then v[landv] = 0 
    437          res=sqrt(u^2+v^2) 
    438          if NOT keyword_set(key_periodic) OR nx NE jpi then res[0,*]=!values.f_nan 
    439          res[*,0]=!values.f_nan 
    440          mask = where(mask eq 0) 
    441          IF mask[0] NE -1 THEN res[mask] = valmask 
     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 
    442444; All kind of average 
    443          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme 
    444          if keyword_set(direc) then res = moyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef) 
    445       END 
    446 ;---------------------------------------------------------------------------- 
    447    endcase 
    448 ;------------------------------------------------------------ 
    449    if keyword_set(key_performance) THEN print, 'temps norme', systime(1)-tempsun 
    450    return, res 
     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 
    451453end 
Note: See TracChangeset for help on using the changeset viewer.