Changeset 167


Ignore:
Timestamp:
09/05/06 14:24:07 (18 years ago)
Author:
smasson
Message:

clean div, curl, grad + minors bugfix

Location:
trunk/SRC
Files:
1 added
4 edited
3 moved

Legend:

Unmodified
Added
Removed

  • trunk/SRC/Computation/curl.pro

    r164 r167  
    1111;  
    1212; @param UU 
    13 ; Matrix representing coordinates of a field of vectors 
     13; Matrix representing the zonal coordinates (U point) of a field of vectors 
     14; A 2D (xy), 3D (xyz or yt) or a structure readable by litchamp and containing 
     15; a 2D (xy), 3D (xyz or yt) array (4d case is not coded yet). 
     16; note that the dimension of the arry must suit the domain dimension. 
    1417; 
    1518; @param VV  
    16 ; Matrix representing coordinates of a field of vectors 
     19; Matrix representing the meridional coordinates (V point) of a field of vectors 
     20; A 2D (xy), 3D (xyz or yt) or a structure readable by litchamp and containing 
     21; a 2D (xy), 3D (xyz or yt) array (4d case is not coded yet). 
     22; note that the dimension of the arry must suit the domain dimension. 
     23; 
     24; @keyword DIREC {type=scalar string} 
     25; Use if you want to call moyenne or grossemoyenne after the div computation 
     26; with a mean done in the DIREC direction 
    1727; 
    1828; @returns RES 
    19 ; A 2d matrix 
     29; the vertical component of the curl of the input data (with the same size) 
    2030; 
    2131; @uses 
    22 ; common.pro 
     32; cm_4cal, cm_4data, cm_4mmesh  
    2333; 
    2434; @restrictions 
    25 ; U and V matrices can be 2 or 4d. 
    26 ; Beware, to discern different configuration of U and V (xy, xyz, xyt, xyzt),  
    27 ; we look at the variable of the common  
    28 ;        -time which contain the calendar in IDL julian days to which U and  
    29 ; V refered to, in the same way as the variable  
    30 ;        -jpt which is the number of time's step to consider in time. 
    31 ; U and V arrays ae cut in the same geographic domain. Because of the gap of  
    32 ; T, U, V and F grids, it is possible that these two arrays has not the same  
    33 ; size and refered to different indexes. In this case, arrays are re-cut on  
    34 ; common indexes and the domain is redefined to match with these common  
    35 ; indexes. To avoid these re-cuts, use the keyword /memeindice in domdef.pro 
    36 ; 
    37 ; 
    38 ; Points on the drawing edge are at !values.f_nan  
     35; 
     36; - Works only for Arakawa C-grid.  
     37; - UU must be on U grid, VV must be on V grid 
     38; - 4d case is not coded yet 
     39; - the common variable jpt is used to differ xyz (jpt=1) and xyt (jpt\=1) cases. 
     40; - U and V arrays are cut in the same geographic domain. Because of the shift between  
     41;   T, U, V and F grids, it is possible that these two arrays do not have the same  
     42;   size and refer to different indexes. In this case, arrays are re-cut on  
     43;   common indexes. To avoid these re-cuts, use the keyword /memeindice in domdef.pro 
     44; - When computing the divergence, we update, vargrid, varname, varunits and the 
     45;   grid position parameters (firstxf, lastxf, nxf, firstyf, lastyf, nyf). 
     46; - points that cannot be computed (domain bondaries, coastline) are set to NaN 
     47; 
     48; @examples 
     49; IDL> @tst_initorca2 
     50; IDL> plt, curl(dist(jpi,jpj), dist(jpi,jpj)) 
    3951; 
    4052; @history  
    4153; Guillaume Roullet (grlod\@ipsl.jussieu.fr) 
    42 ;  
    4354; Sebastien Masson (smasson\@lodyc.jussieu.fr) 
    4455; adaptation to work with a reduce domain 
    … …  
    4859; $Id$ 
    4960; 
     61; @todo  
     62; code the 4d case 
    5063;- 
    5164;------------------------------------------------------------ 
    5265;------------------------------------------------------------ 
    5366;------------------------------------------------------------ 
    54 FUNCTION curl, uu, vv 
     67FUNCTION curl, uu, vv, DIREC = DIREC 
    5568; 
    5669  compile_opt idl2, strictarrsubs 
    5770; 
    58 @common 
    59    tempsun = systime(1)         ; To key_performance 
    60 ; 
    61    IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 
     71@cm_4cal                        ; for jpt 
     72@cm_4data                       ; for varname, vargrid, vardate, varunit, valmask 
     73@cm_4mesh  
     74; 
     75  tempsun = systime(1)          ; To key_performance 
     76; 
     77  IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 
    6278     return, report(['This version of curl is based on Arakawa C-grid.' $ 
    6379                     , 'U and V grids must therefore be defined']) 
    6480; 
    65    u = litchamp(uu) 
    66    v = litchamp(vv) 
    67  
    68    date1 = time[0] 
    69    if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1] 
    70    if (size(u))[0] NE (size(v))[0] then return,  -1 
     81  u = litchamp(uu) 
     82  v = litchamp(vv) 
     83 
     84  szu = size(u) 
     85  szv = size(v) 
     86   
     87  if szu[0] NE szv[0] then return, report('U and V input data must have the same number of dimensions') 
    7188 
    7289;------------------------------------------------------------ 
    7390; We find common points between U and V 
    7491;------------------------------------------------------------ 
    75    indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1] 
    76    indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1] 
    77    indicex = inter(indicexu, indicexv) 
    78    indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1] 
    79    indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1] 
    80    indicey = inter(indiceyu, indiceyv) 
    81    nx = n_elements(indicex)  
    82    ny = n_elements(indicey) 
    83    case 1 of 
    84 ;---------------------------------------------------------------------------- 
     92  indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1] 
     93  indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1] 
     94  indicex = inter(indicexu, indicexv) 
     95  indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1] 
     96  indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1] 
     97  indicey = inter(indiceyu, indiceyv) 
     98  nx = n_elements(indicex)  
     99  ny = n_elements(indicey) 
     100  indice2d = lindgen(jpi, jpj) 
     101  indice2d = indice2d[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1] 
     102;---------------------------------------------------------------------------- 
     103  vargrid = 'F' 
     104  varname = 'vorticity' 
     105  varunits = 's-1' 
     106  if n_elements(valmask) EQ 0 THEN valmask = 1e20 
     107  firstxt = indicex[0] & lastxt = indicex[0]+nx-1 & nxt = nx  
     108  firstyt = indicey[0] & lastyt = indicey[0]+ny-1 & nyt = ny 
     109;---------------------------------------------------------------------------- 
     110;---------------------------------------------------------------------------- 
     111  case 1 of 
    85112;---------------------------------------------------------------------------- 
    86113;xyz 
    87114;---------------------------------------------------------------------------- 
    88       (size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN  
    89          indice2d = lindgen(jpi, jpj) 
    90          indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1] 
     115    szu[0] EQ 3 AND jpt EQ 1:BEGIN  
    91116;------------------------------------------------------------ 
    92117; extraction of U and V on the appropriated domain 
    93118;------------------------------------------------------------ 
    94          case 1 of 
    95             (size(u))[0] NE 3 OR (size(v))[0] NE 3: return,  -1 
    96             (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 
    97              (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
    98             case 1 of 
    99                nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
    100                nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
    101                nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
    102                nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
    103                ELSE : 
    104             endcase 
    105             END 
    106             (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 
    107              (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN  
    108                u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 
    109                v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 
    110             END 
    111             ELSE:return,  -1 
    112          endcase 
    113 ;------------------------------------------------------------ 
    114 ; calculation of the curl 
    115 ;------------------------------------------------------------ 
    116          coefu = (e1u[indice2d])[*]#replicate(1, nzt) 
    117          coefu = reform(coefu, nx, ny, nzt, /over) 
    118          coefu = coefu*(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
    119          terreu = where(coefu EQ 0) 
    120          if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan 
    121           
    122          coefv = (e2v[indice2d])[*]#replicate(1, nzt) 
    123          coefv = reform(coefv, nx, ny, nzt, /over) 
    124          coefv = coefv*(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
    125          terrev = where(coefv EQ 0) 
    126          if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan 
    127  
    128          tabf = (fmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
    129          div = (e1f[indice2d]*e2f[indice2d])[*]#replicate(1, nzt) 
    130          div = reform(div, nx, ny, nzt, /over) 
    131          tabf = tabf/div 
    132 ; 
    133          zu = u*temporary(coefu) 
    134          zv = v*temporary(coefv) 
    135  
    136          psi =  (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0)) 
    137          psi = tabf*psi 
     119      case 1 of 
     120        szu[1] EQ nxu AND szu[2] EQ nyu AND $ 
     121           szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN  
     122          case 1 of 
     123            nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
     124            nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
     125            nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
     126            nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
     127            ELSE : 
     128          endcase 
     129        END 
     130        szu[1] EQ jpi AND szu[2] EQ jpj AND $ 
     131           szv[1] EQ jpi AND szv[2] EQ jpj:BEGIN  
     132          u = u[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *] 
     133          v = v[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *] 
     134        END 
     135        ELSE:return,  -1 
     136      endcase 
     137;------------------------------------------------------------ 
     138; curl computation 
     139;------------------------------------------------------------ 
     140      coefu = ((e1u[indice2d])[*]#replicate(1., nzt)) $ 
     141              * (umask())[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
     142      landu = where(coefu EQ 0) 
     143      if landu[0] NE -1 then coefu[temporary(landu)] = !values.f_nan 
     144       
     145      coefv = ((e2v[indice2d])[*]#replicate(1., nzt)) $ 
     146              *(vmask())[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
     147      landv = where(coefv EQ 0) 
     148      if landv[0] NE -1 then coefv[temporary(landv)] = !values.f_nan 
     149 
     150      tabf = (fmask())[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt] $ 
     151             / ((e1f[indice2d]*e2f[indice2d])[*]#replicate(1., nzt)) 
     152      landf =  where(tabf EQ 0) 
     153; 
     154      zu = temporary(u) * temporary(coefu) 
     155      zv = temporary(v) * temporary(coefv) 
     156 
     157      psi = (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0)) 
     158      psi = temporary(tabf) * temporary(psi) 
    138159;------------------------------------------------------------ 
    139160; Edging put at !values.f_nan 
    140161;------------------------------------------------------------ 
    141          if NOT keyword_set(key_periodic)  OR nx NE jpi then begin 
    142             psi[0, *, *] = !values.f_nan 
    143             psi[nx-1, *, *] = !values.f_nan 
    144          endif 
    145          psi[*, 0, *] = !values.f_nan 
    146          psi[*, ny-1, *] = !values.f_nan 
    147 ; 
    148          if n_elements(valmask) EQ 0 THEN valmask = 1e20 
    149          terref =  where(tabf EQ 0) 
    150          if terref[0] NE -1 then psi[temporary(terref)] = valmask 
    151 ;------------------------------------------------------------ 
    152 ; For the graphic drawing 
    153 ;------------------------------------------------------------ 
    154          domdef, (glagmt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f'] 
    155          if keyword_set(direc) then psi = moyenne(psi,direc,/nan) 
    156  
    157 ;---------------------------------------------------------------------------- 
    158       END 
     162      if NOT keyword_set(key_periodic)  OR nx NE jpi then begin 
     163        psi[0, *, *] = !values.f_nan 
     164        psi[nx-1, *, *] = !values.f_nan 
     165      endif 
     166      psi[*, 0, *] = !values.f_nan 
     167      psi[*, ny-1, *] = !values.f_nan 
     168; 
     169      if landf[0] NE -1 then psi[temporary(landf)] = valmask 
     170      if keyword_set(direc) then psi = moyenne(psi, direc, /nan) 
     171    END 
    159172;---------------------------------------------------------------------------- 
    160173;---------------------------------------------------------------------------- 
    … …  
    162175;---------------------------------------------------------------------------- 
    163176;---------------------------------------------------------------------------- 
    164       date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN  
    165          indice2d = lindgen(jpi, jpj) 
    166          indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1] 
     177    szu[0] EQ 3 AND jpt GT 1:BEGIN  
    167178;------------------------------------------------------------ 
    168179; extraction of U and V on the appropriated domain 
    169180;------------------------------------------------------------ 
    170          case 1 of 
    171             (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 
    172              (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
    173                if nxu NE nx then $ 
    174                 if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]  
    175                IF nxv NE nx THEN $ 
    176                 if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
    177                IF nyu NE ny THEN $ 
    178                 if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]  
    179                IF nyv NE ny THEN $ 
    180                 if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
    181             END 
    182             (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 
    183              (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN  
    184                u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 
    185                v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 
    186             END 
    187             ELSE:BEGIN  
    188                print, 'problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs' 
    189                return, -1 
    190             end 
    191          endcase 
    192 ;---------------------------------------------------------------------------- 
    193 ; Calculation of the curl 
    194 ;---------------------------------------------------------------------------- 
    195          coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt] 
    196          terreu = where(coefu EQ 0) 
    197          if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan 
    198          coefu = temporary(coefu[*])#replicate(1, jpt) 
    199          coefu = reform(coefu, nx, ny, jpt, /over) 
    200 ; 
    201          coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt] 
    202          terrev = where(coefv EQ 0) 
    203          if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan 
    204          coefv = temporary(coefv[*])#replicate(1, jpt) 
    205          coefv = reform(coefv, nx, ny, jpt, /over) 
    206 ; 
    207          tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d]) 
    208          tabf = temporary(tabf[*])#replicate(1, jpt) 
    209          tabf = reform(tabf, nx, ny, jpt, /over) 
    210 ;------------------------------------------------------------ 
    211 ; Calculation of the curl 
    212 ;------------------------------------------------------------ 
    213          zu = u*temporary(coefu) 
    214          zv = v*temporary(coefv) 
    215 ; 
    216          psi =  (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0)) 
    217          psi = tabf*psi 
     181      case 1 of 
     182        szu[1] EQ nxu AND szu[2] EQ nyu AND $ 
     183           szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN  
     184          if nxu NE nx then $ 
     185             if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]  
     186          IF nxv NE nx THEN $ 
     187             if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
     188          IF nyu NE ny THEN $ 
     189             if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]  
     190          IF nyv NE ny THEN $ 
     191             if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
     192        END 
     193        szu[1] EQ jpi AND szu[2] EQ jpj AND $ 
     194           szv[1] EQ jpi AND szv[2] EQ jpj:BEGIN  
     195          u = u[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *] 
     196          v = v[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *] 
     197        END 
     198        ELSE:BEGIN  
     199          print, 'problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs' 
     200          return, -1 
     201        end 
     202      endcase 
     203;---------------------------------------------------------------------------- 
     204; curl computation 
     205;---------------------------------------------------------------------------- 
     206      coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt] 
     207      landu = where(coefu EQ 0) 
     208      if landu[0] NE -1 then coefu[temporary(landu)] = !values.f_nan 
     209      coefu = temporary(coefu[*])#replicate(1., jpt) 
     210; 
     211      coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt] 
     212      landv = where(coefv EQ 0) 
     213      if landv[0] NE -1 then coefv[temporary(landv)] = !values.f_nan 
     214      coefv = temporary(coefv[*])#replicate(1., jpt) 
     215; 
     216      tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d]) 
     217      tabf = temporary(tabf[*])#replicate(1., jpt) 
     218      landf = where(tabf EQ 0) 
     219; 
     220      zu = temporary(u) * temporary(coefu) 
     221      zv = temporary(v) * temporary(coefv) 
     222; 
     223      psi = (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0)) 
     224      psi = temporary(tabf) * temporary(psi) 
    218225;------------------------------------------------------------ 
    219226; extraction of U and V on the appropriated domain 
    220227;------------------------------------------------------------ 
    221          if NOT keyword_set(key_periodic) OR nx NE jpi then begin 
    222             psi[0, *, *] = !values.f_nan 
    223             psi[nx-1, *, *] = !values.f_nan 
    224          endif 
    225          psi[*, 0, *] = !values.f_nan 
    226          psi[*, ny-1, *] = !values.f_nan 
    227          if n_elements(valmask) EQ 0 THEN valmask = 1e20 
    228          terref =  where(tabf EQ 0) 
    229          if terref[0] NE -1 then psi[temporary(terref)] = valmask 
    230 ;---------------------------------------------------------------------------- 
    231          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f'] 
    232          if keyword_set(direc) then psi = grossemoyenne(psi,direc,/nan) 
    233 ;---------------------------------------------------------------------------- 
    234       END 
     228      if NOT keyword_set(key_periodic) OR nx NE jpi then begin 
     229        psi[0, *, *] = !values.f_nan 
     230        psi[nx-1, *, *] = !values.f_nan 
     231      endif 
     232      psi[*, 0, *] = !values.f_nan 
     233      psi[*, ny-1, *] = !values.f_nan 
     234      if landf[0] NE -1 then psi[temporary(landf)] = valmask 
     235      if keyword_set(direc) then psi = grossemoyenne(psi, direc, /nan) 
     236    END 
    235237;---------------------------------------------------------------------------- 
    236238;---------------------------------------------------------------------------- 
    … …  
    238240;---------------------------------------------------------------------------- 
    239241;---------------------------------------------------------------------------- 
    240       date1 NE date2 AND (size(u))[0] EQ 4:BEGIN  
    241          return, report('non code!') 
    242       END 
     242    szu[0] EQ 4:BEGIN  
     243      return, report('Case not coded contact saxo team or make a do loop!') 
     244    END 
    243245;---------------------------------------------------------------------------- 
    244246;---------------------------------------------------------------------------- 
    … …  
    246248;---------------------------------------------------------------------------- 
    247249;---------------------------------------------------------------------------- 
    248       ELSE:BEGIN                ;xy 
    249          indice2d = lindgen(jpi, jpj) 
    250          indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1] 
    251 ;------------------------------------------------------------ 
    252 ;------------------------------------------------------------ 
    253          case 1 of 
    254             (size(u))[0] NE 2 OR (size(v))[0] NE 2: return,  -1 
    255             (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 
    256              (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
    257                if nxu NE nx then $ 
    258                 if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]  
    259                IF nxv NE nx THEN $ 
    260                 if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 
    261                IF nyu NE ny THEN $ 
    262                 if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]  
    263                IF nyv NE ny THEN $ 
    264                 if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 
    265             END 
    266             (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 
    267              (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN  
    268                u = u[indice2d] 
    269                v = v[indice2d] 
    270             END 
    271             ELSE:return,  -1 
    272          endcase 
    273 ;------------------------------------------------------------ 
    274 ; Calculation of the curl 
    275 ;------------------------------------------------------------ 
    276          coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt] 
    277          terreu = where(coefu EQ 0) 
    278          if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan 
    279          coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt] 
    280          terrev = where(coefv EQ 0) 
    281          if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan 
    282          tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d]) 
    283 ; 
    284          zu = u*temporary(coefu) 
    285          zv = v*temporary(coefv) 
    286  
    287          psi =  (shift(zv, -1, 0)-zv) + (zu-shift(zu, 0, -1)) 
    288          psi = tabf*psi 
    289  
     250    szu[0] EQ 2:BEGIN  
     251;------------------------------------------------------------ 
     252;------------------------------------------------------------ 
     253      case 1 of 
     254        szu[1] EQ nxu AND szu[2] EQ nyu AND $ 
     255           szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN  
     256          if nxu NE nx then $ 
     257             if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]  
     258          IF nxv NE nx THEN $ 
     259             if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 
     260          IF nyu NE ny THEN $ 
     261             if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]  
     262          IF nyv NE ny THEN $ 
     263             if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 
     264        END 
     265        szu[1] EQ jpi AND szu[2] EQ jpj AND $ 
     266           szv[1] EQ jpi AND szv[2] EQ jpj:BEGIN  
     267          u = u[indice2d] 
     268          v = v[indice2d] 
     269        END 
     270        ELSE:return,  -1 
     271      endcase 
     272;------------------------------------------------------------ 
     273; curl computation 
     274;------------------------------------------------------------ 
     275      coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt] 
     276      landu = where(coefu EQ 0) 
     277      if landu[0] NE -1 then coefu[temporary(landu)] = !values.f_nan 
     278      coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt] 
     279      landv = where(coefv EQ 0) 
     280      if landv[0] NE -1 then coefv[temporary(landv)] = !values.f_nan 
     281      tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d]) 
     282      landf =  where(tabf EQ 0) 
     283; 
     284      zu = temporary(u) * temporary(coefu) 
     285      zv = temporary(v) * temporary(coefv) 
     286 
     287      psi = (shift(zv, -1, 0)-zv) + (zu-shift(zu, 0, -1)) 
     288      psi = temporary(tabf) * temporary(psi) 
    290289;------------------------------------------------------------ 
    291290; Edging put at !values.f_nan 
    292291;------------------------------------------------------------ 
    293          if  NOT keyword_set(key_periodic) OR nx NE jpi then begin 
    294             psi[0, *] = !values.f_nan 
    295             psi[nx-1, *] = !values.f_nan 
    296          endif 
    297          psi[*, 0] = !values.f_nan 
    298          psi[*, ny-1] = !values.f_nan 
    299 ; 
    300          if n_elements(valmask) EQ 0 THEN valmask = 1e20 
    301          terref =  where(tabf EQ 0) 
    302          if terref[0] NE -1 then psi[temporary(terref)] = valmask 
    303 ;------------------------------------------------------------ 
    304 ; for the graphic drawing 
    305 ;------------------------------------------------------------ 
    306          domdef, (glamt[indice2d])[0, 0], (glamf[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphif[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f'] 
    307          if keyword_set(direc) then psi = moyenne(psi,direc,/nan) 
    308  
    309       END 
    310 ;---------------------------------------------------------------------------- 
    311    endcase 
    312 ;------------------------------------------------------------ 
    313    if keyword_set(key_performance) THEN print, 'temps curl', systime(1)-tempsun  
    314 ;------------------------------------------------------------ 
    315    vargrid = 'F' 
    316    varname = 'vorticity' 
    317    return, psi 
     292      if  NOT keyword_set(key_periodic) OR nx NE jpi then begin 
     293        psi[0, *] = !values.f_nan 
     294        psi[nx-1, *] = !values.f_nan 
     295      endif 
     296      psi[*, 0] = !values.f_nan 
     297      psi[*, ny-1] = !values.f_nan 
     298; 
     299      if landf[0] NE -1 then psi[temporary(landf)] = valmask 
     300      if keyword_set(direc) then psi = moyenne(psi, direc, /nan) 
     301    END 
     302;---------------------------------------------------------------------------- 
     303;---------------------------------------------------------------------------- 
     304    ELSE:return, report('U and V input arrays must have 2, 3 or 4 dimensions') 
     305  ENDCASE 
     306;------------------------------------------------------------ 
     307  if keyword_set(key_performance) THEN print, 'temps curl', systime(1)-tempsun  
     308 
     309  return, psi 
    318310end 

  • trunk/SRC/Computation/div.pro

    r164 r167  
    55; 
    66; @file_comments 
    7 ; calculation of the divergence of a 2d field 
     7; compute the horizontal divergence of a vectors field 
    88; 
    99; @categories 
    … …  
    1111; 
    1212; @param UU  
    13 ; Matrix representing coordinates of a field of vectors 
     13; Matrix representing the zonal coordinates (U point) of a field of vectors 
     14; A 2D (xy), 3D (xyz or yt) or a structure readable by litchamp and containing 
     15; a 2D (xy), 3D (xyz or yt) array (4d case is not coded yet). 
     16; note that the dimension of the arry must suit the domain dimension. 
    1417; 
    1518; @param VV  
    16 ; Matrix representing coordinates of a field of vectors 
     19; Matrix representing the meridional coordinates (V point) of a field of vectors 
     20; A 2D (xy), 3D (xyz or yt) or a structure readable by litchamp and containing 
     21; a 2D (xy), 3D (xyz or yt) array (4d case is not coded yet). 
     22; note that the dimension of the arry must suit the domain dimension. 
     23; 
     24; @keyword DIREC {type=scalar string} 
     25; Use if you want to call moyenne or grossemoyenne after the div computation 
     26; (stupid ?) with a mean done in the DIREC direction 
    1727; 
    1828; @returns RES 
    19 ; A 2d matrix 
     29; the divergence of the input data (with the same size) 
    2030; 
    2131; @uses 
    22 ; common.pro 
     32; cm_4cal, cm_4data, cm_4mmesh  
    2333; 
    2434; @restrictions 
    25 ; U and V matrices can be 2 or 4d. 
    26 ; Beware, to discern different configuration of U and V (xy, xyz, xyt, xyzt),  
    27 ; we look at the variable of the common  
    28 ;        -time which contain the calendar in IDL julian days to which U and  
    29 ; V refered to, in the same way as the variable  
    30 ;        -jpt which is the number of time's step to consider in time. 
    31 ; U and V arrays are cut in the same geographic domain. Because of the gap of  
    32 ; T, U, V and F grids, it is possible that these two arrays has not the same  
    33 ; size and refered to different indexes. In this case, arrays are re-cut on  
    34 ; common indexes and the domain is redefined to match with these common  
    35 ; indexes. To avoid these re-cuts, use the keyword /memeindice in domdef.pro 
    36 ; 
    37 ; 
    38 ; Points on the drawing edge are at !values.f_nan  
    39 ; 
     35; 
     36; - Works only for Arakawa C-grid.  
     37; - UU must be on U grid, VV must be on V grid 
     38; - 4d case is not coded yet 
     39; - the common variable jpt is used to differ xyz (jpt=1) and xyt (jpt\=1) cases. 
     40; - U and V arrays are cut in the same geographic domain. Because of the shift between  
     41;   T, U, V and F grids, it is possible that these two arrays do not have the same  
     42;   size and refer to different indexes. In this case, arrays are re-cut on  
     43;   common indexes. To avoid these re-cuts, use the keyword /memeindice in domdef.pro 
     44; - When computing the divergence, we update, vargrid, varname, varunits and the 
     45;   grid position parameters (firstxt, lastxt, nxt, firstyt, lastyt, nyt). 
     46; - points that cannot be computed (domain bondaries, coastline) are set to NaN 
     47; 
     48; @examples 
     49; IDL> @tst_initorca2 
     50; IDL> plt, div(dist(jpi,jpj), dist(jpi,jpj)) 
    4051; 
    4152; @history 
    42 ; Guillaume Roullet (grlod\@ipsl.jussieu.fr) 
    43 ;                      Creation : printemps 1998 
    44 ;                      Sebastien Masson (smasson\@lodyc.jussieu.fr) 
    45 ;                      adaptation to work with a reduce domain 
    46 ;                      12/1/2000; 
     53; Guillaume Roullet (grlod\@ipsl.jussieu.fr): creation; spring 1998 
     54; Sebastien Masson (smasson\@lodyc.jussieu.fr) 
     55; adaptation to work with a reduce domain; 12/1/2000 
    4756; 
    4857; @version 
    4958; $Id$ 
    5059; 
     60; @todo  
     61; code the 4d case 
    5162;- 
    5263;------------------------------------------------------------ 
    5364;------------------------------------------------------------ 
    5465;------------------------------------------------------------ 
    55 FUNCTION div, uu, vv 
     66FUNCTION div, uu, vv, DIREC = DIREC 
    5667; 
    5768  compile_opt idl2, strictarrsubs 
    5869; 
    59    tempsun = systime(1)         ; For key_performance 
    60 @common 
    61 ; 
    62    IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 
     70@cm_4cal                        ; for jpt 
     71@cm_4data                       ; for varname, vargrid, vardate, varunit, valmask 
     72@cm_4mesh  
     73; 
     74  tempsun = systime(1)          ; For key_performance 
     75; 
     76  IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 
    6377     return, report(['This version of div is based on Arakawa C-grid.' $ 
    6478                     , 'U and V grids must therefore be defined']) 
    6579; 
    66    u = litchamp(uu) 
    67    v = litchamp(vv) 
    68 ; 
    69    date1 = time[0] 
    70    if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1] 
    71    if (size(u))[0] NE (size(v))[0] then return,  -1 
     80  u = litchamp(uu) 
     81  v = litchamp(vv) 
     82; 
     83  szu = size(u) 
     84  szv = size(v) 
     85   
     86  if szu[0] NE szv[0] then return, report('U and V input data must have the same number of dimensions') 
    7287 
    7388;------------------------------------------------------------ 
    7489; We find common points between U and V 
    7590;------------------------------------------------------------ 
    76    indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1] 
    77    indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1] 
    78    indicex = inter(indicexu, indicexv) 
    79    indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1] 
    80    indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1] 
    81    indicey = inter(indiceyu, indiceyv) 
    82    nx = n_elements(indicex)  
    83    ny = n_elements(indicey) 
    84    indice2d = lindgen(jpi, jpj) 
    85    indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1] 
    86    case 1 of 
    87 ;---------------------------------------------------------------------------- 
     91  indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1] 
     92  indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1] 
     93  indicex = inter(indicexu, indicexv) 
     94  indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1] 
     95  indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1] 
     96  indicey = inter(indiceyu, indiceyv) 
     97  nx = n_elements(indicex)  
     98  ny = n_elements(indicey) 
     99  indice2d = lindgen(jpi, jpj) 
     100  indice2d = indice2d[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1] 
     101;---------------------------------------------------------------------------- 
     102  vargrid = 'T' 
     103  varname = 'div' 
     104  varunits = '1.e6*s-1' 
     105  if n_elements(valmask) EQ 0 THEN valmask = 1.e20 
     106  firstxt = indicex[0] & lastxt = indicex[0]+nx-1 & nxt = nx  
     107  firstyt = indicey[0] & lastyt = indicey[0]+ny-1 & nyt = ny 
     108;---------------------------------------------------------------------------- 
     109;---------------------------------------------------------------------------- 
     110  case 1 of 
    88111;---------------------------------------------------------------------------- 
    89112;xyz 
    90113;---------------------------------------------------------------------------- 
    91       (size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN  
     114    szu[0] EQ 3 AND jpt EQ 1:BEGIN  
    92115;------------------------------------------------------------ 
    93116; extraction of U and V on the appropriated domain 
    94117;------------------------------------------------------------ 
    95          case 1 of 
    96             (size(v))[0] NE 3: return,  -1 
    97             (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 
    98              (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
    99                case 1 of 
    100                   nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
    101                   nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
    102                   nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
    103                   nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
    104                   ELSE : 
    105                endcase 
    106             END 
    107             (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 
    108              (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN  
    109                u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 
    110                v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 
    111             END 
    112             ELSE:BEGIN  
    113                zdiv = -1 
    114                GOTO, sortie 
    115             end 
    116              
    117          endcase 
    118 ;------------------------------------------------------------ 
    119 ; calcul de la divergence 
    120 ;------------------------------------------------------------ 
    121          zu = (e2u[indice2d])[*]#replicate(1, nzt) 
    122          zu = reform(zu, nx, ny, nzt, /over) 
    123          zu = temporary(u)*temporary(zu) $ 
    124           *(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
    125          terreu = where(zu EQ 0) 
    126          if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan 
    127 ; 
    128          zv = (e1v[indice2d])[*]#replicate(1, nzt) 
    129          zv = reform(zv, nx, ny, nzt, /over) 
    130          zv = temporary(v)*temporary(zv) $ 
    131           *(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
    132          terrev = where(zv EQ 0) 
    133          if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan 
    134 ; 
    135          zdiv = 1e6/(e1t[indice2d]*e2t[indice2d]) 
    136          zdiv = (zdiv)[*]#replicate(1, nzt) 
    137          zdiv = reform(zdiv, nx, ny, nzt, /over) 
    138          zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0)) $ 
    139           *tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
     118      case 1 of 
     119        szu[1] EQ nxu AND szu[2] EQ nyu AND $ 
     120           szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN  
     121          case 1 of 
     122            nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
     123            nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
     124            nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
     125            nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
     126            ELSE : 
     127          endcase 
     128        END 
     129        szu[1] EQ jpi AND szu[2] EQ jpj AND $ 
     130           szv[1] EQ jpi AND szv[2] EQ jpj:BEGIN  
     131          u = u[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *] 
     132          v = v[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *] 
     133        END 
     134        ELSE:return, -1 
     135      endcase 
     136;------------------------------------------------------------ 
     137; divergence computation 
     138;------------------------------------------------------------ 
     139      zu = (e2u[indice2d])[*]#replicate(1., nzt) 
     140      landu = where((umask())[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0) 
     141      if landu[0] NE -1 then zu[temporary(landu)] = !values.f_nan 
     142      zu = temporary(u) * temporary(zu) 
     143; 
     144      zv = (e1v[indice2d])[*]#replicate(1., nzt) 
     145      landv = where((vmask())[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0) 
     146      if landv[0] NE -1 then zv[temporary(landv)] = !values.f_nan 
     147      zv = temporary(v) * temporary(zv) 
     148; 
     149      zdiv = (e1t[indice2d]*e2t[indice2d])[*]#replicate(1.e6, nzt) 
     150      zdiv = ( zu - shift(zu, 1, 0, 0) + zv - shift(zv, 0, 1, 0) ) * temporary(zdiv) 
    140151;------------------------------------------------------------ 
    141152; Edging put at !values.f_nan 
    142153;------------------------------------------------------------ 
    143          if  NOT keyword_set(key_periodic) OR nx NE jpi then begin 
    144             zdiv[0, *, *] = !values.f_nan 
    145             zdiv[nx-1, *, *] = !values.f_nan 
    146          endif 
    147          zdiv[*, 0, *] = !values.f_nan 
    148          zdiv[*, ny-1, *] = !values.f_nan 
    149 ; 
    150          zdiv = temporary(zdiv) 
    151 ; 
    152          if n_elements(valmask) EQ 0 THEN valmask = 1e20 
    153          terre =  where(tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0) 
    154          if terre[0] NE -1 then zdiv[temporary(terre)] = valmask 
    155 ;------------------------------------------------------------ 
    156 ; For the graphic drawing 
    157 ;------------------------------------------------------------ 
    158          vargrid = 'T' 
    159          varname = 'div' 
    160          varunits = '1e6*s-1' 
    161          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t'] 
    162          if keyword_set(direc) then  zdiv = moyenne(zdiv,direc,/nan) 
    163       END 
     154      if  NOT keyword_set(key_periodic) OR nx NE jpi then begin 
     155        zdiv[0, *, *] = !values.f_nan 
     156        zdiv[nx-1, *, *] = !values.f_nan 
     157      endif 
     158      zdiv[*, 0, *] = !values.f_nan 
     159      zdiv[*, ny-1, *] = !values.f_nan 
     160; 
     161      land = where(tmask[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0) 
     162      if land[0] NE -1 then zdiv[temporary(land)] = valmask 
     163      if keyword_set(direc) then  zdiv = moyenne(zdiv, direc, /nan) 
     164    END 
    164165;---------------------------------------------------------------------------- 
    165166;---------------------------------------------------------------------------- 
    … …  
    167168;---------------------------------------------------------------------------- 
    168169;---------------------------------------------------------------------------- 
    169       date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN  
     170    szu[0] EQ 3 AND jpt GT 1:BEGIN  
    170171;------------------------------------------------------------ 
    171172; extraction of U and V on the appropriated domain 
    172173;------------------------------------------------------------ 
    173          case 1 of 
    174             (size(u))[0] NE 3 OR (size(v))[0] NE 3: return,  -1 
    175             (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 
    176              (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
    177                case 1 of 
    178                   nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
    179                   nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
    180                   nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
    181                   nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
    182                   ELSE : 
    183                endcase 
    184             END 
    185             (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 
    186              (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN  
    187                u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 
    188                v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 
    189             END 
    190             ELSE:return,  -1 
    191          endcase 
    192 ;------------------------------------------------------------ 
    193 ; Calculation of the divergence 
    194 ;------------------------------------------------------------ 
    195          zu = e2u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt] 
    196          terreu = where(zu EQ 0) 
    197          if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan 
    198          zu = (zu)[*]#replicate(1, jpt) 
    199          zu = reform(zu, nx, ny, jpt, /over) 
    200          zu = temporary(u)*temporary(zu) 
    201 ; 
    202          zv = e1v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt] 
    203          terrev = where(zv EQ 0) 
    204          if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan 
    205          zv = (zv)[*]#replicate(1, jpt) 
    206          zv = reform(zv, nx, ny, jpt, /over) 
    207          zv = temporary(v)*temporary(zv) 
    208 ; 
    209          zdiv = 1e6*tmask[indice2d+jpi*jpj*firstzt]/(e1t[indice2d]*e2t[indice2d]) 
    210          zdiv = (zdiv)[*]#replicate(1, jpt) 
    211          zdiv = reform(zdiv, nx, ny, jpt, /over) 
    212          terre =  where(zdiv EQ 0) 
    213          zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0)) 
     174      case 1 of 
     175        szu[1] EQ nxu AND szu[2] EQ nyu AND $ 
     176           szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN  
     177          case 1 of 
     178            nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
     179            nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
     180            nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
     181            nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
     182            ELSE : 
     183          endcase 
     184        END 
     185        szu[1] EQ jpi AND szu[2] EQ jpj AND $ 
     186           szv[1] EQ jpi AND szv[2] EQ jpj:BEGIN  
     187          u = u[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *] 
     188          v = v[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *] 
     189        END 
     190        ELSE:return, -1 
     191      endcase 
     192;------------------------------------------------------------ 
     193; divergence computation 
     194;------------------------------------------------------------ 
     195      zu = e2u[indice2d] 
     196      landu = where((umask())[indice2d+jpi*jpj*firstzt] EQ 0) 
     197      if landu[0] NE -1 then zu[temporary(landu)] = !values.f_nan 
     198      zu = (temporary(zu))[*]#replicate(1., jpt) 
     199      zu = temporary(u) * temporary(zu) 
     200; 
     201      zv = e1v[indice2d] 
     202      landv = where((vmask())[indice2d+jpi*jpj*firstzt] EQ 0) 
     203      if landv[0] NE -1 then zv[temporary(landv)] = !values.f_nan 
     204      zv = (temporary(zv))[*]#replicate(1., jpt) 
     205      zv = temporary(v) * temporary(zv) 
     206; 
     207      zdiv = (e1t[indice2d]*e2t[indice2d])[*]#replicate(1.e6, jpt) 
     208      zdiv = ( zu - shift(zu, 1, 0, 0) + zv - shift(zv, 0, 1, 0) ) * temporary(zdiv) 
    214209;------------------------------------------------------------ 
    215210; Edging put at !values.f_nan 
    216211;------------------------------------------------------------ 
    217          if  NOT keyword_set(key_periodic) OR nx NE jpi then begin 
    218             zdiv[0, *, *] = !values.f_nan 
    219             zdiv[nx-1, *, *] = !values.f_nan 
    220          endif 
    221          zdiv[*, 0, *] = !values.f_nan 
    222          zdiv[*, ny-1, *] = !values.f_nan 
    223 ; 
    224          if n_elements(valmask) EQ 0 THEN valmask = 1e20 
    225          if terre[0] NE -1 then zdiv[temporary(terre)] = valmask 
    226 ;------------------------------------------------------------ 
    227 ; for the graphic drawing 
    228 ;------------------------------------------------------------ 
    229          vargrid = 'T' 
    230          varname = 'div' 
    231          varunits = '1e6*s-1' 
    232          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t'] 
    233          if keyword_set(direc) then  zdiv = grossemoyenne(zdiv,direc,/nan) 
    234       END 
     212      if  NOT keyword_set(key_periodic) OR nx NE jpi then begin 
     213        zdiv[0, *, *] = !values.f_nan 
     214        zdiv[nx-1, *, *] = !values.f_nan 
     215      endif 
     216      zdiv[*, 0, *] = !values.f_nan 
     217      zdiv[*, ny-1, *] = !values.f_nan 
     218; 
     219      land = where(tmask[indice2d+jpi*jpj*firstzt] EQ 0, cnt) 
     220      if land[0] NE -1 then BEGIN  
     221        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*lindgen(jpt)) 
     222        zdiv[temporary(land)] = valmask 
     223      ENDIF 
     224      if keyword_set(direc) then  zdiv = grossemoyenne(zdiv, direc, /nan) 
     225    END 
    235226;---------------------------------------------------------------------------- 
    236227;---------------------------------------------------------------------------- 
    … …  
    238229;---------------------------------------------------------------------------- 
    239230;---------------------------------------------------------------------------- 
    240       date1 NE date2 AND (size(u))[0] EQ 4:BEGIN  
    241          return, report('non code!') 
    242       END 
     231    szu[0] EQ 4:BEGIN  
     232      return, report('Case not coded contact saxo team or make a do loop!') 
     233    END 
    243234;---------------------------------------------------------------------------- 
    244235;---------------------------------------------------------------------------- 
    … …  
    246237;---------------------------------------------------------------------------- 
    247238;---------------------------------------------------------------------------- 
    248       ELSE:BEGIN                ;xy 
    249          indice3d = lindgen(jpi, jpj, jpk) 
    250          indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt] 
     239    szu[0] EQ 2:BEGIN  
    251240;------------------------------------------------------------ 
    252241; extraction of U and V on the appropriated domain 
    253242;------------------------------------------------------------ 
    254          case 1 of 
    255             (size(u))[0] NE 2 OR (size(v))[0] NE 2: BEGIN  
    256                zdiv = -1 
    257                GOTO, sortie 
    258             end 
    259             (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 
    260              (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
    261                case 1 of 
    262                   nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 
    263                   nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 
    264                   nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 
    265                   nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 
    266                   ELSE : 
    267                endcase 
    268             END 
    269             (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 
    270              (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN  
    271                u = u[indice2d] 
    272                v = v[indice2d] 
    273             END 
    274             ELSE:return,  -1 
    275          endcase 
    276 ;------------------------------------------------------------ 
    277 ; Calculation of the divergence 
    278 ;------------------------------------------------------------ 
    279          zu = temporary(u)*e2u[indice2d]*(umask())[indice3d] 
    280          terreu = where(zu EQ 0) 
    281          if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan 
    282          zv = temporary(v)*e1v[indice2d]*(vmask())[indice3d] 
    283          terrev = where(zv EQ 0) 
    284          if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan 
    285          zdiv = zu-shift(zu, 1, 0)+zv-shift(zv, 0, 1) 
    286          zdiv = temporary(zdiv)*tmask[indice3d]/(e1t[indice2d]*e2t[indice2d]) 
     243      case 1 of 
     244        szu[1] EQ nxu AND szu[2] EQ nyu AND $ 
     245           szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN  
     246          case 1 of 
     247            nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 
     248            nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 
     249            nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 
     250            nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 
     251            ELSE : 
     252          endcase 
     253        END 
     254        szu[1] EQ jpi AND szu[2] EQ jpj AND $ 
     255           szv[1] EQ jpi AND szv[2] EQ jpj:BEGIN  
     256          u = u[indice2d] 
     257          v = v[indice2d] 
     258        END 
     259        ELSE:return, -1 
     260      endcase 
     261;------------------------------------------------------------ 
     262; divergence computation 
     263;------------------------------------------------------------ 
     264      zu = e2u[indice2d] 
     265      landu = where((umask())[indice2d+jpi*jpj*firstzt] EQ 0) 
     266      if landu[0] NE -1 then zu[temporary(landu)] = !values.f_nan 
     267      zu = temporary(u) * temporary(zu) 
     268 
     269      zv = e1v[indice2d] 
     270      landv = where((vmask())[indice2d+jpi*jpj*firstzt] EQ 0) 
     271      if landv[0] NE -1 then zv[temporary(landv)] = !values.f_nan 
     272      zv = temporary(v) * temporary(zv) 
     273 
     274      zdiv = 1.e6 / (e1t[indice2d]*e2t[indice2d]) 
     275      zdiv = ( zu - shift(zu, 1, 0) + zv - shift(zv, 0, 1) ) * temporary(zdiv) 
    287276;------------------------------------------------------------ 
    288277; Edging put at !values.f_nan 
    289278;------------------------------------------------------------ 
    290          if  NOT keyword_set(key_periodic) OR nx NE jpi then begin 
    291             zdiv[0, *] = !values.f_nan 
    292             zdiv[nx-1, *] = !values.f_nan 
    293          endif 
    294          zdiv[*, 0] = !values.f_nan 
    295          zdiv[*, ny-1] = !values.f_nan 
    296 ; 
    297          zdiv = temporary(zdiv)*1e6 
    298 ; 
    299          if n_elements(valmask) EQ 0 THEN valmask = 1e20 
    300          terre =  where(tmask[indice3d] EQ 0) 
    301          if terre[0] NE -1 then zdiv[temporary(terre)] = valmask 
    302 ;------------------------------------------------------------ 
    303 ; for the graphic drawing 
    304 ;------------------------------------------------------------ 
    305          vargrid = 'T' 
    306          varname = 'div' 
    307          varunits = '1e6*s-1' 
    308          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t'] 
    309          if keyword_set(direc) then  zdiv = moyenne(zdiv,direc,/nan) 
    310  
    311       END 
    312    endcase 
    313     
    314 ;------------------------------------------------------------ 
    315 sortie: 
    316    if keyword_set(key_performance) THEN print, 'temps div', systime(1)-tempsun  
    317     
    318    return, zdiv 
     279      if  NOT keyword_set(key_periodic) OR nx NE jpi then begin 
     280        zdiv[0, *] = !values.f_nan 
     281        zdiv[nx-1, *] = !values.f_nan 
     282      endif 
     283      zdiv[*, 0] = !values.f_nan 
     284      zdiv[*, ny-1] = !values.f_nan 
     285; 
     286      land =  where(tmask[indice2d+jpi*jpj*firstzt] EQ 0) 
     287      if land[0] NE -1 then zdiv[temporary(land)] = valmask 
     288      if keyword_set(direc) then zdiv = moyenne(zdiv, direc, /nan) 
     289    END 
     290;---------------------------------------------------------------------------- 
     291;---------------------------------------------------------------------------- 
     292    ELSE:return, report('U and V input arrays must have 2, 3 or 4 dimensions') 
     293  ENDCASE    
     294;------------------------------------------------------------ 
     295  if keyword_set(key_performance) THEN print, 'temps div', systime(1)-tempsun  
     296   
     297  return, zdiv 
    319298end 

  • trunk/SRC/Computation/grad.pro

    r164 r167  
    44;+ 
    55; @file_comments 
    6 ; 
     6; compute the gradient of a variable 
    77; 
    88; @categories 
    9 ; 
     9; Calculation 
    1010; 
    1111; @param FIELD 
    12 ; 
    13 ; 
    14 ; @param DIREC 
    15 ; 
    16 ; 
    17 ; @returns 
    18 ; 
     12; The field for which we want to compute the gradient.  A 2D (xy), 
     13; 3D (xyz or yt) or 4D (xyzt) array or a structure readable by litchamp 
     14; and containing a 2D (xy), 3D (xyz or yt) or 4D (xyzt) array. 
     15; note that the dimension of the arry must suit the domain dimension. 
     16; 
     17; @param DIREC {type=scalar string} 
     18; the gradient direction: 'x', 'y', 'z' 
     19; 
     20; @returns RES {type=2D, 3D or 4D array} 
     21; the gradient of the input data (with the same size) 
    1922; 
    2023; @uses 
    21 ; 
     24; cm_4cal, cm_4data, cm_4mmesh  
    2225; 
    2326; @restrictions 
    24 ; 
     27; - Works only for Arakawa C-grid.  
     28; - When computing the gradient, the result is not on the same grid point 
     29;   than the input data. In consequence, we update, vargrid and the grid position 
     30;   parameters (firstx[tuvf], lastx[tuvf], nx[tuvf], firsty[tuvf], lasty[tuvf], 
     31;   ny[tuvf], firstz[tw], lastz[tw], nz[tw]). 
     32; - points that cannot be computed (domain bondaries, coastline) are set to NaN 
     33; - the common variable jpt is used to differ xyz (jpt=1) and xyt (jpt\=1) cases. 
    2534; 
    2635; @examples 
    27 ; 
     36; IDL> @tst_initorca2 
     37; IDL> plt, grad({arr:gphit,g:'T'}, 'x') 
     38; IDL> plt, grad({arr:gphit,g:'T'}, 'y') 
    2839; 
    2940; @history 
    … …  
    3344; $Id$ 
    3445; 
    35 ; @todo  
    36 ; seb: remplir les truc! 
    3746;- 
    3847;------------------------------------------------------------ 
    … …  
    4352  compile_opt idl2, strictarrsubs 
    4453; 
    45 @common 
    46 ;------------------------------------------------------------ 
    47 ; 
    48    IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 
     54@cm_4cal   ; for jpt 
     55@cm_4data  ; for varname, vargrid, vardate, varunit, valmask 
     56@cm_4mesh  
     57;------------------------------------------------------------ 
     58; 
     59  IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 
    4960     return, report(['This version of grad is based on Arakawa C-grid.' $ 
    5061                     , 'U and V grids must therefore be defined']) 
    5162; 
    52    res = litchamp(field) 
    53    taille = size(res) 
    54    grille, mask, glam, gphi, gdep, nx, ny, nz $ 
    55            , firstx, firsty, firstz, lastx, lasty, lastz      
    56    if n_elements(valmask) EQ 0 then valmask = 1e20 
    57    case strupcase(vargrid) of 
    58       'T':BEGIN 
    59          case direc of 
    60             'x':BEGIN  
    61                divi = e1u[firstx:lastx, firsty:lasty] 
    62                newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz] 
    63                vargrid = 'U' 
    64                firstxu = firstxt & lastxu = lastxt & nxu = nxt  
    65                firstyu = firstyt & lastyu = lastyt & nyu = nyt 
    66             END 
    67             'y':BEGIN 
    68                divi = e2v[firstx:lastx, firsty:lasty] 
    69                newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz] 
    70                vargrid = 'V' 
    71                firstxv = firstxt & lastxv = lastxt & nxv = nxt  
    72                firstyv = firstyt & lastyv = lastyt & nyv = nyt 
    73             END 
    74             'z':BEGIN 
    75                divi = e3w[firstz:lastz] 
    76                newmask = mask 
    77                vargrid = 'W' 
    78                firstzw = firstzt & lastzw = lastzt & nzw = nzt  
    79             END 
    80             ELSE:return, report('Bad definition of direction argument') 
    81          ENDCASE 
    82       END 
    83       'W':BEGIN 
    84          case direc of 
    85             'x':BEGIN  
    86                divi = e1u[firstx:lastx, firsty:lasty] 
    87                newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz] 
    88                vargrid = 'U' 
    89                firstxu = firstxt & lastxu = lastxt & nxu = nxt  
    90                firstyu = firstyt & lastyu = lastyt & nyu = nyt 
    91             END 
    92             'y':BEGIN 
    93                divi = e2v[firstx:lastx, firsty:lasty] 
    94                newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz] 
    95                vargrid = 'V' 
    96                firstxv = firstxt & lastxv = lastxt & nxv = nxt  
    97                firstyv = firstyt & lastyv = lastyt & nyv = nyt 
    98             END 
    99             'z':BEGIN 
    100                divi = e3t[firstz:lastz] 
    101                newmask = mask 
    102                vargrid = 'T' 
    103                firstzt = firstzw & lastzt = lastzw & nzt = nzw  
    104             END 
    105             ELSE:return, report('Bad definition of direction argument') 
    106          endcase 
    107       END 
    108       'U':BEGIN 
    109          case direc of 
    110             'x':BEGIN 
    111                divi = (shift(e1t, -1, 0))[firstx:lastx, firsty:lasty] 
    112                newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz] 
    113                vargrid = 'T' 
    114                firstxt = firstxu & lastxt = lastxu & nxt = nxu  
    115                firstyt = firstyu & lastyt = lastyu & nyt = nyu 
    116             END 
    117             'y':BEGIN 
    118                divi = e2f[firstx:lastx, firsty:lasty] 
    119                newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz] 
    120                vargrid = 'F' 
    121                firstxf = firstxu & lastxf = lastxu & nxf = nxu  
    122                firstyf = firstyu & lastyf = lastyu & nyf = nyu 
    123             END 
    124             'z':BEGIN 
    125                divi = e3w[firstz:lastz] 
    126                newmask = mask 
    127                vargrid = 'W' 
    128                firstzw = firstzt & lastzw = lastzt & nzw = nzt  
    129             END 
    130             ELSE:return, report('Bad definition of direction argument') 
    131          endcase 
    132       END 
    133       'V':BEGIN 
    134          case direc of 
    135             'x':BEGIN 
    136                divi = e1f[firstx:lastx, firsty:lasty] 
    137                newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz] 
    138                vargrid = 'F' 
    139                firstxf = firstxv & lastxf = lastxv & nxf = nxv 
    140                firstyf = firstyv & lastyf = lastyv & nyf = nyv 
    141             END 
    142             'y':BEGIN 
    143                divi = (shift(e2t, 0, -1))[firstx:lastx, firsty:lasty] 
    144                newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz] 
    145                vargrid = 'T' 
    146                firstxt = firstxv & lastxt = lastxv & nxt = nxv 
    147                firstyt = firstyv & lastyt = lastyv & nyt = nyv 
    148             END 
    149             'z':BEGIN 
    150                divi = e3w[firstz:lastz] 
    151                newmask = mask 
    152                vargrid = 'W' 
    153                firstzw = firstzt & lastzw = lastzt & nzw = nzt  
    154             END 
    155             ELSE:return, report('Bad definition of direction argument') 
    156          endcase 
    157       END 
    158 ;       'F':BEGIN 
     63  res = litchamp(field) 
     64  szres = size(res) 
     65  grille, mask, glam, gphi, gdep, nx, ny, nz $ 
     66          , firstx, firsty, firstz, lastx, lasty, lastz      
     67; 
     68  if n_elements(valmask) EQ 0 then valmask = 1.e20 
     69  varname = 'grad of '+varname 
     70  varunit = varunit+'/m' 
     71  case strupcase(vargrid) of 
     72    'T':BEGIN 
     73      case direc of 
     74        'x':BEGIN  
     75          divi = e1u[firstx:lastx, firsty:lasty] 
     76          newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz] 
     77          vargrid = 'U' 
     78          firstxu = firstxt & lastxu = lastxt & nxu = nxt  
     79          firstyu = firstyt & lastyu = lastyt & nyu = nyt 
     80        END 
     81        'y':BEGIN 
     82          divi = e2v[firstx:lastx, firsty:lasty] 
     83          newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz] 
     84          vargrid = 'V' 
     85          firstxv = firstxt & lastxv = lastxt & nxv = nxt  
     86          firstyv = firstyt & lastyv = lastyt & nyv = nyt 
     87        END 
     88        'z':BEGIN 
     89          divi = e3w[firstz:lastz] 
     90          newmask = mask 
     91          vargrid = 'W' 
     92          firstzw = firstzt & lastzw = lastzt & nzw = nzt  
     93        END 
     94        ELSE:return, report('Bad definition of direction argument') 
     95      ENDCASE 
     96    END 
     97    'W':BEGIN 
     98      case direc of 
     99        'x':BEGIN  
     100          divi = e1u[firstx:lastx, firsty:lasty] 
     101          newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz] 
     102          vargrid = 'U' 
     103          firstxu = firstxt & lastxu = lastxt & nxu = nxt  
     104          firstyu = firstyt & lastyu = lastyt & nyu = nyt 
     105        END 
     106        'y':BEGIN 
     107          divi = e2v[firstx:lastx, firsty:lasty] 
     108          newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz] 
     109          vargrid = 'V' 
     110          firstxv = firstxt & lastxv = lastxt & nxv = nxt  
     111          firstyv = firstyt & lastyv = lastyt & nyv = nyt 
     112        END 
     113        'z':BEGIN 
     114          divi = e3t[firstz:lastz] 
     115          newmask = mask 
     116          vargrid = 'T' 
     117          firstzt = firstzw & lastzt = lastzw & nzt = nzw  
     118        END 
     119        ELSE:return, report('Bad definition of direction argument') 
     120      endcase 
     121    END 
     122    'U':BEGIN 
     123      case direc of 
     124        'x':BEGIN 
     125          divi = (shift(e1t, -1, 0))[firstx:lastx, firsty:lasty] 
     126          newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz] 
     127          vargrid = 'T' 
     128          firstxt = firstxu & lastxt = lastxu & nxt = nxu  
     129          firstyt = firstyu & lastyt = lastyu & nyt = nyu 
     130        END 
     131        'y':BEGIN 
     132          divi = e2f[firstx:lastx, firsty:lasty] 
     133          newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz] 
     134          vargrid = 'F' 
     135          firstxf = firstxu & lastxf = lastxu & nxf = nxu  
     136          firstyf = firstyu & lastyf = lastyu & nyf = nyu 
     137        END 
     138        'z':BEGIN 
     139          divi = e3w[firstz:lastz] 
     140          newmask = mask 
     141          vargrid = 'W' 
     142          firstzw = firstzt & lastzw = lastzt & nzw = nzt  
     143        END 
     144        ELSE:return, report('Bad definition of direction argument') 
     145      endcase 
     146    END 
     147    'V':BEGIN 
     148      case direc of 
     149        'x':BEGIN 
     150          divi = e1f[firstx:lastx, firsty:lasty] 
     151          newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz] 
     152          vargrid = 'F' 
     153          firstxf = firstxv & lastxf = lastxv & nxf = nxv 
     154          firstyf = firstyv & lastyf = lastyv & nyf = nyv 
     155        END 
     156        'y':BEGIN 
     157          divi = (shift(e2t, 0, -1))[firstx:lastx, firsty:lasty] 
     158          newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz] 
     159          vargrid = 'T' 
     160          firstxt = firstxv & lastxt = lastxv & nxt = nxv 
     161          firstyt = firstyv & lastyt = lastyv & nyt = nyv 
     162        END 
     163        'z':BEGIN 
     164          divi = e3w[firstz:lastz] 
     165          newmask = mask 
     166          vargrid = 'W' 
     167          firstzw = firstzt & lastzw = lastzt & nzw = nzt  
     168        END 
     169        ELSE:return, report('Bad definition of direction argument') 
     170      endcase 
     171    END 
     172    'F':BEGIN 
    159173;          case direc of 
    160174;             'x':divi = (shift(e1v, -1, 0))[firstx:lastx, firsty:lasty] 
    … …  
    163177;             ELSE:return, report('Bad definition of direction argument') 
    164178;          endcase 
    165 ;       END 
    166       ELSE:return, report('Bad definition of vargrid') 
    167    ENDCASE 
    168    res = fitintobox(res) 
    169    IF n_elements(res) EQ 1 AND res[0] EQ -1 THEN return, res 
    170    case 1 of 
     179      return, report('F grid: case not coded, please contact SAXO team...') 
     180    END 
     181    ELSE:return, report('Bad definition of vargrid') 
     182  ENDCASE 
     183  res = fitintobox(temporary(res)) 
     184  IF n_elements(res) EQ 1 AND res[0] EQ -1 THEN return, res 
     185  case 1 of 
    171186;---------------------------------------------------------------------------- 
    172187;xy 
    173188;---------------------------------------------------------------------------- 
    174       taille[0] EQ 2:BEGIN 
    175          earth = where(mask[*, *, firstz] EQ 0) 
    176          if earth[0] NE -1 then res[earth] = !values.f_nan 
    177          case direc of 
    178             'x':BEGIN  
    179                res = (shift(res, -1, 0)-res)/divi 
    180                if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *] = !values.f_nan 
    181                if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0) 
    182             END 
    183             'y':BEGIN  
    184                res = (shift(res, 0, -1)-res)/divi 
    185                res[*, ny-1] = !values.f_nan 
    186                if vargrid EQ 'T' OR vargrid EQ 'U' then res =   shift(res, 0, 1)            
    187             END 
    188             ELSE:return,  report('Bad definition of direction argument for the type of array') 
    189          ENDCASE 
    190          earth = where(newmask[*, *, firstz] EQ 0) 
    191          if earth[0] NE -1 then res[earth] = valmask 
    192       END 
     189    szres[0] EQ 2:BEGIN 
     190      land = where((temporary(mask))[*, *, firstz] EQ 0) 
     191      if land[0] NE -1 then res[temporary(land)] = !values.f_nan 
     192      case direc of 
     193        'x':BEGIN  
     194          res = (shift(res, -1, 0)-res)/temporary(divi) 
     195          if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *] = !values.f_nan 
     196          if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0) 
     197        END 
     198        'y':BEGIN  
     199          res = (shift(res, 0, -1)-res)/temporary(divi) 
     200          res[*, ny-1] = !values.f_nan 
     201          if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1)            
     202        END 
     203        ELSE:return,  report('Bad definition of direction argument for the type of array') 
     204      ENDCASE 
     205      land = where((temporary(newmask))[*, *, firstz] EQ 0) 
     206      if land[0] NE -1 then res[temporary(land)] = valmask 
     207    END 
    193208;---------------------------------------------------------------------------- 
    194209;xyt 
    195210;---------------------------------------------------------------------------- 
    196       taille[0] EQ 3 AND jpt NE 1:BEGIN 
    197          earth = where(mask[*, *, firstz] EQ 0) 
    198          if earth[0] NE -1 then BEGIN 
    199             earth = earth#replicate(1, jpt)+replicate(1, n_elements(earth))#(nx*ny*lindgen(jpt)) 
    200             res[earth[*]] = !values.f_nan 
    201          ENDIF 
    202          divi = divi[*]#replicate(1, jpt) 
    203          case direc of 
    204             'x':BEGIN  
    205                res = (shift(res, -1, 0, 0)-res)/divi 
    206                if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan 
    207                if vargrid EQ 'T' OR vargrid EQ 'V'  then res = shift(res, 1, 0, 0) 
    208             END 
    209             'y':BEGIN  
    210                res = (shift(res, 0, -1, 0)-res)/divi 
    211                res[*, ny-1, *] = !values.f_nan 
    212                if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0)         
    213             END 
    214             ELSE:return,  report('Bad definition of direction argument for the type of array') 
    215          ENDCASE 
    216          earth = where(newmask[*, *, firstz] EQ 0) 
    217          if earth[0] NE -1 THEN res[earth] = valmask 
    218       END 
     211    szres[0] EQ 3 AND jpt NE 1:BEGIN 
     212      land = where((temporary(mask))[*, *, firstz] EQ 0, cnt) 
     213      if land[0] NE -1 then BEGIN 
     214        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*lindgen(jpt)) 
     215        res[temporary(land)] = !values.f_nan 
     216      ENDIF 
     217      divi = (temporary(divi))[*]#replicate(1., jpt) 
     218      case direc of 
     219        'x':BEGIN  
     220          res = (shift(res, -1, 0, 0)-res)/temporary(divi) 
     221          if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan 
     222          if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0) 
     223        END 
     224        'y':BEGIN  
     225          res = (shift(res, 0, -1, 0)-res)/temporary(divi) 
     226          res[*, ny-1, *] = !values.f_nan 
     227          if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0)         
     228        END 
     229        ELSE:return,  report('Bad definition of direction argument for the type of array') 
     230      ENDCASE 
     231      land = where((temporary(newmask))[*, *, firstz] EQ 0, cnt) 
     232      if land[0] NE -1 then BEGIN 
     233        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*lindgen(jpt)) 
     234        res[temporary(land)] = valmask 
     235      ENDIF 
     236    END 
    219237;---------------------------------------------------------------------------- 
    220238;xyz 
    221239;---------------------------------------------------------------------------- 
    222       taille[0] EQ 3 AND jpt EQ 1:BEGIN 
    223          earth = where(mask EQ 0) 
    224          if earth[0] NE -1 then res[earth] = !values.f_nan 
    225          case direc OF 
    226             'x':BEGIN  
    227                divi = divi[*]#replicate(1, nz) 
    228                res = (shift(res, -1, 0, 0)-res)/divi 
    229                if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan 
    230                if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0) 
    231             END 
    232             'y':BEGIN  
    233                divi = divi[*]#replicate(1, nz) 
    234                res = (shift(res, 0, -1, 0)-res)/divi 
    235                res[*, ny-1, *] = !values.f_nan 
    236                if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0)         
    237             END 
    238             'z':BEGIN 
    239                divi = reform(replicate(1, nx*ny)#divi, nx, ny, nz) 
    240                if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz) 
    241                if vargrid EQ 'W' THEN BEGIN  
    242                   res = (shift(res, 0, 0, 1)-res)/divi 
    243                   res[*, *, 0] = !values.f_nan 
    244                ENDIF ELSE BEGIN 
    245                   res = (res-shift(res, 0, 0, -1))/divi 
    246                   res[*, *, nz-1] = !values.f_nan 
    247                ENDELSE 
    248                if earth[0] NE -1 then res[earth] = valmask 
    249             END 
    250          ENDCASE 
    251       END 
    252 ;------------------------------------------------------------ 
     240    szres[0] EQ 3 AND jpt EQ 1:BEGIN 
     241      land = where(mask EQ 0) 
     242      if land[0] NE -1 then res[temporary(land)] = !values.f_nan 
     243      case direc OF 
     244        'x':BEGIN  
     245          divi = (temporary(divi))[*]#replicate(1., nz) 
     246          res = (shift(res, -1, 0, 0)-res)/temporary(divi) 
     247          if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan 
     248          if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0) 
     249        END 
     250        'y':BEGIN  
     251          divi = (temporary(divi))[*]#replicate(1., nz) 
     252          res = (shift(res, 0, -1, 0)-res)/temporary(divi) 
     253          res[*, ny-1, *] = !values.f_nan 
     254          if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0)         
     255        END 
     256        'z':BEGIN 
     257          divi = replicate(1., nx*ny)#(temporary(divi))[*] 
     258          if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, /overwrite) 
     259          if vargrid EQ 'W' THEN BEGIN  
     260            res = (shift(res, 0, 0, 1)-res)/temporary(divi) 
     261            res[*, *, 0] = !values.f_nan 
     262          ENDIF ELSE BEGIN 
     263            res = (res-shift(res, 0, 0, -1))/temporary(divi) 
     264            res[*, *, nz-1] = !values.f_nan 
     265          ENDELSE 
     266        END 
     267      ENDCASE 
     268      land = where(temporary(newmask) EQ 0) 
     269      if land[0] NE -1 then res[temporary(land)] = valmask 
     270    END 
    253271;---------------------------------------------------------------------------- 
    254272;xyzt 
    255273;---------------------------------------------------------------------------- 
    256       taille[0] EQ 4:BEGIN 
    257          earth = where((mask[*]#replicate(1, jpt)) EQ 0) 
    258          if earth[0] NE -1 then res[earth] = !values.f_nan 
    259          case direc OF 
    260             'x':BEGIN  
    261                divi = divi[*]#replicate(1, nz*jpt) 
    262                res = (shift(res, -1, 0, 0, 0)-res)/divi 
    263                if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *, *] = !values.f_nan 
    264                if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0, 0) 
    265             END 
    266             'y':BEGIN  
    267                divi = divi[*]#replicate(1, nz*jpt) 
    268                res = (shift(res, 0, -1, 0, 0)-res)/divi 
    269                res[*, ny-1, *, *] = !values.f_nan 
    270                if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0, 0)         
    271             END 
    272             'z':BEGIN 
    273                divi = replicate(1, nx*ny)#divi 
    274                divi = reform(divi[*]#replicate(1, jpt), nx, ny, nz, jpt, /over) 
    275                if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, jpt) 
    276                if vargrid EQ 'W' THEN BEGIN  
    277                   res = (shift(res, 0, 0, 1, 0)-res)/divi 
    278                   res[*, *, 0, *] = !values.f_nan 
    279                ENDIF ELSE BEGIN 
    280                   res = (res-shift(res, 0, 0, -1, 0))/divi 
    281                   res[*, *, nz-1, *] = !values.f_nan 
    282                ENDELSE 
    283             END 
    284          ENDCASE 
    285          if earth[0] NE -1 then res[earth] = valmask 
    286       END 
    287 ;------------------------------------------------------------ 
    288 ;------------------------------------------------------------ 
    289    endcase 
    290    varname = 'grad of '+varname 
    291    varunit = varunit+'/m' 
    292  
    293  
    294 ;------------------------------------------------------------ 
    295    return, res 
    296 end 
    297  
    298  
    299  
    300  
    301  
    302  
     274    szres[0] EQ 4:BEGIN 
     275      land = where(temporary(mask) EQ 0, cnt) 
     276      if land[0] NE -1 then BEGIN 
     277        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*nz*lindgen(jpt)) 
     278        res[temporary(land)] = !values.f_nan 
     279      ENDIF 
     280      case direc OF 
     281        'x':BEGIN  
     282          divi = (temporary(divi))[*]#replicate(1., nz*jpt) 
     283          res = (shift(res, -1, 0, 0, 0)-res)/temporary(divi) 
     284          if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *, *] = !values.f_nan 
     285          if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0, 0) 
     286        END 
     287        'y':BEGIN  
     288          divi = (temporary(divi))[*]#replicate(1., nz*jpt) 
     289          res = (shift(res, 0, -1, 0, 0)-res)/temporary(divi) 
     290          res[*, ny-1, *, *] = !values.f_nan 
     291          if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0, 0)         
     292        END 
     293        'z':BEGIN 
     294          divi = replicate(1., nx*ny)#(temporary(divi))[*] 
     295          divi = (temporary(divi))[*]#replicate(1L, jpt) 
     296          if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, jpt, /overwrite) 
     297          if vargrid EQ 'W' THEN BEGIN  
     298            res = (shift(res, 0, 0, 1, 0)-res)/temporary(divi) 
     299            res[*, *, 0, *] = !values.f_nan 
     300          ENDIF ELSE BEGIN 
     301            res = (res-shift(res, 0, 0, -1, 0))/temporary(divi) 
     302            res[*, *, nz-1, *] = !values.f_nan 
     303          ENDELSE 
     304        END 
     305      ENDCASE 
     306      land = where(newmask EQ 0, cnt) 
     307      if land[0] NE -1 then BEGIN 
     308        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*nz*lindgen(jpt)) 
     309        res[temporary(land)] = valmask 
     310      ENDIF 
     311    END 
     312;------------------------------------------------------------ 
     313;------------------------------------------------------------ 
     314    ELSE:return, report('input array must have 2, 3 or 4 dimensions') 
     315  ENDCASE 
     316;------------------------------------------------------------ 
     317 
     318 
     319;------------------------------------------------------------ 
     320  return, res 
     321END 
     322 
     323 
     324 
     325 
     326 
     327 

  • trunk/SRC/Postscript/openps.pro

    r136 r167  
    128128           , LANDSCAPE = 1 - key_portrait, PORTRAIT = key_portrait $ 
    129129           , xsize = xs, ysize = ys, xoffset = xoff, yoffset = yoff $ 
    130            , bits_per_pixel = 8, _extra = ex 
     130           , bits_per_pixel = 8, language_level = 2, _extra = ex 
    131131; to make smaller postcripts 
    132132   IF NOT (keyword_set(keeppfont) OR keyword_set(keep_pfont)) $ 

  • trunk/SRC/ToBeReviewed/INIT/initncdf.pro

    r163 r167  
    123123  if keyword_set(xaxisname) then xaxisname = strlowcase(xaxisname) ELSE xaxisname = 'x' 
    124124  xvarid = where(namevar EQ xaxisname OR namevar EQ 'longitude' $ 
    125                  OR namevar EQ 'nav_lon' OR namevar EQ 'lon') 
     125                 OR namevar EQ 'nav_lon' OR namevar EQ 'lon' $ 
     126                 OR namevar EQ 'NbLongitudes') 
    126127  xvarid = xvarid[0] 
    127128  if xvarid EQ -1 then begin 
    … …  
    145146; find the yaxis 
    146147  if keyword_set(yaxisname) then yaxisname = strlowcase(yaxisname) ELSE yaxisname = 'y' 
    147   yvarid = where(namevar EQ yaxisname OR namevar EQ 'latitude' OR namevar EQ 'nav_lat' OR namevar EQ 'lat') 
     148  yvarid = where(namevar EQ yaxisname OR namevar EQ 'latitude' $ 
     149                 OR namevar EQ 'nav_lat' OR namevar EQ 'lat' $ 
     150                 OR namevar EQ 'NbLatitudes') 
    148151  yvarid = yvarid[0] 
    149152  if yvarid EQ -1 then begin 

  • trunk/SRC/ToBeReviewed/WIDGET/AUTOUR_de_XXX/scanfile.pro

    r150 r167  
    114114    namedim[dimiq] = strlowcase(tmpname) 
    115115  ENDFOR 
    116 ; we are looking for a x dimension... 
    117   dimidx = where(namedim EQ 'x' OR strmid(namedim, 0, 3) EQ 'lon' OR strmid(namedim, 0, 3) EQ 'xi_' OR namedim EQ 'xt_i7_156') 
    118   dimidx = dimidx[0] 
    119   if dimidx EQ -1 then begin 
    120     print, 'one of the dimensions must have the name: ''x'' or ''lon...'' or ''xi_...'' or ''xt_i7_156''' 
    121     stop 
    122   endif 
    123 ; we are looking for a y dimension... 
    124   dimidy = where(namedim EQ 'y' OR strmid(namedim, 0, 3) EQ 'lat' OR strmid(namedim, 4) EQ 'eta_' OR namedim EQ 'yt_j6_75') 
    125   dimidy = dimidy[0] 
    126   if dimidy EQ -1 then begin 
    127     print, 'one of the dimensions must have the name: ''y'' or ''lat...'' or ''eta_...'' or ''yt_j6_75''' 
    128     stop 
    129   endif 
     116; we are looking for a x dimension with a name matching one of the following regular expression: 
     117  testname = ['longitude', 'lon', 'x', 'longitude*', 'lon*', 'x*', '*longitude*', '*lon*', '*x*'] 
     118  cnt = -1 
     119  ii = 0 
     120  WHILE cnt NE 1 AND ii LT n_elements(testname) DO BEGIN 
     121    dimidx = where(strmatch(namedim, testname[ii]) EQ 1, cnt) 
     122    ii = ii+1 
     123  ENDWHILE 
     124  CASE cnt OF 
     125    0:begin 
     126      dummy = report(['none of the dimensions name matches one of the following regular expression:' $ 
     127                      , '''longitude'', ''lon'', ''x'', ''longitude*'', ''lon*'', ''x*'', ''*longitude*'', ''*lon*'', ''*x*''' $ 
     128                      , ' => we cannot find the x dimension']) 
     129      stop 
     130    END 
     131    1:dimidx = dimidx[0] 
     132    ELSE:begin 
     133      dummy = report(['several (and not one unique) dimensions name matches the following regular expression:' $ 
     134                      , '''longitude'', ''lon'', ''x'', ''longitude*'', ''lon*'', ''x*'', ''*longitude*'', ''*lon*'', ''*x*''' $ 
     135                      , ' => we cannot find the x dimension']) 
     136      stop 
     137    END 
     138  ENDCASE 
     139; we are looking for a y dimension with a name matching one of the following regular expression: 
     140  testname = ['latitude', 'lat', 'y', 'latitude*', 'lat*', 'y*', 'eta_*', '*latitude*', '*lat*', '*y*'] 
     141  cnt = -1 
     142  ii = 0 
     143  WHILE cnt NE 1 AND ii LT n_elements(testname) DO BEGIN 
     144    dimidy = where(strmatch(namedim, testname[ii]) EQ 1, cnt) 
     145    ii = ii+1 
     146  ENDWHILE 
     147  CASE cnt OF 
     148    0:begin 
     149      dummy = report(['none of the dimensions name matches one of the following regular expression:' $ 
     150                      , '''latitude'', ''lat'', ''y'', ''latitude*'', ''lat*'', ''y*'', ''eta_*'', ''*latitude*'', ''*lat*'', ''*y*''' $ 
     151                      , ' => we cannot find the y dimension']) 
     152      stop 
     153    END 
     154    1:dimidy = dimidy[0] 
     155    ELSE:begin 
     156      dummy = report(['several (and not one unique) dimensions name matches the following regular expression:' $ 
     157                      , '''latitude'', ''lat'', ''y'', ''latitude*'', ''lat*'', ''y*'', ''eta_*'', ''*latitude*'', ''*lat*'', ''*y*''' $ 
     158                      , ' => we cannot find the x dimension']) 
     159      stop 
     160    END 
     161  ENDCASE 
    130162;------------------------------------------------------------ 
    131163; name of all variables 
    … …  
    156188    varid = 0 
    157189    repeat BEGIN 
    158       invar = ncdf_varinq(cdfid, varid)  
    159       varid = varid+1 
    160     endrep until n_elements(invar.dim) EQ 1 AND invar.dim[0] EQ infile.recdim 
     190      IF varid LT infile.nvars THEN BEGIN 
     191        invar = ncdf_varinq(cdfid, varid)  
     192        varid = varid+1         
     193      ENDIF ELSE varid = 0 
     194    endrep until (n_elements(invar.dim) EQ 1 AND invar.dim[0] EQ infile.recdim) OR (varid EQ 0) 
    161195    varid = varid-1 
    162196; 
    … …  
    179213        for attiq = 0, invar.natts-1 do attnames[attiq] = ncdf_attname(cdfid, varid, attiq) 
    180214        if (where(attnames EQ 'units'))[0] EQ -1 then BEGIN 
    181           dummy = report('Attribut ''units'' not found for the variable '+varid.name+'!C we create a fake calendar ...') 
     215          dummy = report('Attribut ''units'' not found for the variable '+invar.name+'!C we create a fake calendar ...') 
    182216          fakecal = 1 
    183217          time = date0fk + lindgen(jpt) 
    … …  
    218252          ENDELSE 
    219253; 
    220 ; BEWARE we have to recuperate the calendar attribute and ajust TIME by consequence... 
    221 ; 
    222 ; 
    223 ; We pass TIME in IDL julian days 
     254; BEWARE we have to get back the calendar attribute and ajust time by consequence... 
     255; 
     256; 
     257; We pass time in IDL julian days 
    224258; 
    225259          unite = strlowcase(unite) 

  • trunk/SRC/ToBeReviewed/WIDGET/COMPOUND_WIDGET/cw_domain.pro

    r157 r167  
    662662; 
    663663  min2 = gdep2[indice1] 
    664   if indice2 EQ jpk-1 then BEGIN  
    665     max2 = max([gdept, gdepw]) 
    666     max2 = strtrim(string(max2, format = '(e8.0)'), 1) 
    667     max2 = float('1'+strmid(max2, 1))+float(max2) 
    668   ENDIF ELSE max2 = gdep1[indice2+1] 
    669   if max2 EQ min2 then max2 = min2+1 
     664  if indice2 EQ jpk-1 then max2 = ceil(max([gdept, gdepw])) $ 
     665  ELSE max2 = gdep1[indice2+1] 
     666  if floor(max2) LE floor(min2) then max2 = min2+1 
    670667  rien = cw_slider_pm(basez, value = (min2 > boxzoom[4]) > boxzoom[5] < max2 $ 
    671668                      , uvalue = {name:'depth2'}, minimum = min2, maximum =  max2 $ 
Note: See TracChangeset for help on using the changeset viewer.