Ignore:
Timestamp:
04/28/06 14:18:03 (18 years ago)
Author:
pinsard
Message:

upgrade of GRILLE/Utilities according to cerbere.lodyc.jussieu.fr: /usr/home/smasson/SAXO_RD/ : files

Location:
trunk/ToBeReviewed/GRILLE
Files:
13 copied

Legend:

Unmodified
Added
Removed

  • trunk/ToBeReviewed/GRILLE/changegrid.pro

    r12 r13  
    2828   if newgrid.filetype EQ 'batch file' THEN BEGIN  
    2929      createpro, '@'+strmid(newgrid.filename[0], 0, strlen(newgrid.filename)-4) $ 
    30           , filename = isadirectory(io = homedir, title = 'Bad definition of Homedir') $ 
    31           +'for_createpro.pro' 
     30          , filename = myuniquetmpdir +'for_createpro.pro' 
    3231      return, 1 
    3332   ENDIF ELSE BEGIN  
    … …  
    3635; 
    3736; 
    38    key_periodique = newgrid.key_periodique 
     37   key_periodic = newgrid.key_periodic 
     38;    
     39   @updateold 
    3940   domdef 
    40    if newgrid.triangulation EQ 1 then triangles = triangule() ELSE triangles = -1 
     41   if newgrid.triangulation EQ 1 then triangles_list = triangule() $ 
     42   ELSE triangles_list = -1 
     43; 
     44  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     45   @updateold 
     46  ENDIF 
    4147; 
    4248   return, 1 

  • trunk/ToBeReviewed/GRILLE/cmpgrid.pro

    r12 r13  
    1414; we compare the structure which caracterise the grid whith 
    1515; ccmeshparameters 
     16; 
    1617; 
    1718   case 1 of 

  • trunk/ToBeReviewed/GRILLE/domdef.pro

    r12 r13  
    1010; CATEGORY: 
    1111; 
    12 ; CALLING SEQUENCE:domdef [,lon1, lon2, lat1, lat2[,prof1,prof2]] ou 
     12; CALLING SEQUENCE:domdef [,lon1, lon2, lat1, lat2[,vert1,vert2]] ou 
    1313;             bien domdef,vecteur 
    1414; 
    1515; INPUTS:(facultatif), [vecteur a] 2, 4 ou 6 elements:; 
    1616; sans l''utilisation des mots cles index,xindex,yindex,zindex: 
    17 ; *prof1, prof2: pour un domaine 3D dont la partie horizontale couvre tout 
     17; *vert1, vert2: pour un domaine 3D dont la partie horizontale couvre tout 
    1818;  glam et gphi 
    1919; *lon1, lon2, lat1, lat2: 
    2020; definissant les longitudes min. max et les latitudes min, max du  domaine a 
    2121; etudier (tous les niveaux sont selectiones) 
    22 ; *lon1,lon2,lat1,lat2,prof1,prof2 pour specifier les profondeurs. 
     22; *lon1,lon2,lat1,lat2,vert1,vert2 pour specifier les profondeurs. 
    2323; 
    2424; KEYWORD PARAMETERS: 
     25; 
     26;       ENDPOINTS: a four elements vector [x1,y1,x2,y2] used to specify 
     27;       that domdef must define the box used to make a plot (pltz, pltt, 
     28;       plt1d) done strictly along the line (that can have any direction) 
     29;       starting at (x1, y1) ending at (x2, y2). When defining endpoints, 
     30;       you must also define TYPE which define the type of plots 
     31;       ('pltz', 'xt', 'yt', 'zt', 'x', 'y', 'z', 't') will used 
     32;       ENDPOINTS keywords 
    2533; 
    2634;       FINDALWAYS:oblige a redefinir une boite meme qd auqun point 
    … …  
    2836;       grille. 
    2937; 
    30 ;       GRILLE:un string ou un vecteur de string contennant le nom des 
     38;       GRIDTYPE:un string ou un vecteur de string contennant le nom des 
    3139;       grilles (determinees uniquement par : 'T','U','V','W','F') pour 
    3240;       lesquelles le calcul doit etre fait. par 
    … …  
    6169;   -nzt,w:entier qui contient le nombre de pts en profondeur de 
    6270;    la grille reduite au domaine 3D 
    63 ;   -premierxt,u,v,f: le premier indice qui delimite le sous domaine 
     71;   -firstxt,u,v,f: le first indice qui delimite le sous domaine 
    6472;   suivant x 
    65 ;   -premieryt,u,v,f: le premier indice qui delimite le sous domaine 
     73;   -firstyt,u,v,f: le first indice qui delimite le sous domaine 
    6674;   suivant y 
    67 ;   -premierzt,w: le premier indice qui delimite le sous domaine 
     75;   -firstzt,w: le first indice qui delimite le sous domaine 
    6876;   suivant z 
    69 ;   -dernierxt,u,v,f: le dernier indice qui delimite le sous domaine 
     77;   -lastxt,u,v,f: le last indice qui delimite le sous domaine 
    7078;   suivant x 
    71 ;   -dernieryt,u,v,f: le dernier indice qui delimite le sous domaine 
     79;   -lastyt,u,v,f: le last indice qui delimite le sous domaine 
    7280;   suivant y 
    73 ;   -dernierzt,w: le dernier indice qui delimite le sous domaine 
     81;   -lastzt,w: le last indice qui delimite le sous domaine 
    7482;   suivant z 
    7583; 
    … …  
    8593; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr) 
    8694;                       8/2/98 
     95; rewrite everything, debug and spee-up Sebastien Masson April 2005 
    8796;- 
    8897;------------------------------------------------------------ 
    8998;------------------------------------------------------------ 
    9099;------------------------------------------------------------ 
    91 pro domdef,x1,x2,y1,y2,z1,z2, FINDALWAYS = findalways, GRILLE=grille, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex 
    92 @common 
    93    tempsun = systime(1)         ; pour key_performance 
    94 ;------------------------------------------------------------------- 
    95 ; determination en fonction du nombre d''arguments de 
    96 ; lon1,lon2,lat1,lat2,prof1,prof2 
    97 ;------------------------------------------------------------------- 
    98    if not keyword_set(grille) then grille=['T', 'U', 'V', 'W', 'F'] ELSE $  
    99     for i = 0, n_elements(grille)-1 do grille[i] = strupcase(grille[i]) 
    100    if keyword_set(memeindices) then grille=['T', grille] 
    101    CASE N_PARAMS() OF 
    102       0: BEGIN 
     100pro domdef, x1, x2, y1, y2, z1, z2, FINDALWAYS = findalways $ 
     101            , GRIDTYPE = gridtype, MEMEINDICES = memeindices $ 
     102            , XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex $ 
     103            , ENDPOINTS = endpoints, TYPE = type $ 
     104            , INDEX = index, _extra = ex 
     105;------------------------------------------------------------ 
     106; include commons 
     107@cm_4mesh 
     108  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     109@updatenew 
     110@updatekwd 
     111  ENDIF 
     112;--------------------- 
     113  tempsun = systime(1)          ; pour key_performance 
     114; 
     115  CASE N_PARAMS() OF 
     116    0: 
     117    1: 
     118    2: 
     119    4: 
     120    6: 
     121    ELSE:BEGIN 
     122      ras = report('Bad number of parameter in the call of domdef.') 
     123      RETURN 
     124    END 
     125  ENDCASE 
     126; 
     127  IF keyword_set(endpoints) THEN BEGIN 
     128    IF NOT keyword_set(type) THEN BEGIN  
     129      dummy = report(['If domdef is used do find the box associated' $ 
     130                      , 'to endpoints, you must also specify type keyword']) 
     131      return 
     132    ENDIF  
     133    CASE N_PARAMS() OF 
     134      0: 
     135      1:boxzoom = [x1] 
     136      2:boxzoom = [x1, x2] 
     137      4:boxzoom = [x1, x2, y1, y2] 
     138      6:boxzoom = [x1, x2, y1, y2, z1, z2] 
     139    ENDCASE 
     140    section, BOXZOOM = boxzoom, ENDPOINTS = endpoints, TYPE = type, /ONLYBOX 
     141    return 
     142  ENDIF    
     143 
     144;------------------------------------------------------------------- 
     145; recall domdef when there is only one input parameter 
     146;------------------------------------------------------------------- 
     147; 
     148  IF N_PARAMS() EQ 1 THEN BEGIN 
     149    CASE n_elements(x1) OF 
     150      2:domdef, x1[0], x1[1], FINDALWAYS = findalways, GRIDTYPE = gridtype, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex 
     151      4:domdef, x1[0], x1[1], x1[2], x1[3], FINDALWAYS = findalways, GRIDTYPE = gridtype, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex 
     152      6:domdef, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5], FINDALWAYS = findalways, GRIDTYPE = gridtype, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex 
     153      ELSE:BEGIN 
     154        ras = report('Bad number of elements in x1') 
     155        RETURN 
     156      END 
     157    ENDCASE 
     158    RETURN   
     159  ENDIF 
     160;------------------------------------------------------------------- 
     161; default definitions and checks 
     162;------------------------------------------------------------------- 
     163  IF NOT keyword_set(gridtype) THEN gridtype = ['T', 'U', 'V', 'W', 'F'] $ 
     164  ELSE gridtype = strupcase(gridtype) 
     165  IF keyword_set(memeindices) THEN gridtype = ['T', gridtype] 
     166; default definitions 
     167  lon1t = 99999. & lon2t = -99999. & lat1t = 99999. & lat2t = -99999. 
     168  lon1u = 99999. & lon2u = -99999. & lat1u = 99999. & lat2u = -99999. 
     169  lon1v = 99999. & lon2v = -99999. & lat1v = 99999. & lat2v = -99999. 
     170  lon1f = 99999. & lon2f = -99999. & lat1f = 99999. & lat2f = -99999. 
     171  vert1t = 99999. & vert2t = -99999. & vert1w = 99999. & vert2w = -99999. 
     172; 
     173  IF N_PARAMS() EQ 2 THEN GOTO, vertical 
     174; 
     175;------------------------------------------------------------------- 
     176;------------------------------------------------------------------- 
     177; define all horizontal parameters ... 
    103178; lon1 et lon2 
    104          test = 1 
    105          IF (where(grille eq 'T'))[0] NE -1 THEN test = [test, glamt[*]] 
    106          IF (where(grille eq 'W'))[0] NE -1 AND (where(grille eq 'T'))[0] EQ -1 THEN $ 
    107           test = [test, glamt[*]] 
    108          IF (where(grille eq 'U'))[0] NE -1 AND n_elements(glamu) NE 0 THEN test = [test, glamu[*]] 
    109          IF (where(grille eq 'V'))[0] NE -1 AND n_elements(glamv) NE 0 THEN test = [test, glamv[*]] 
    110          IF (where(grille eq 'F'))[0] NE -1 AND n_elements(glamf) NE 0 THEN test = [test, glamf[*]] 
    111          test = test[1:n_elements(test)-1]  
    112          lon1=min([test], max = max) 
    113          lon2=max 
    114179; lat1 et lat2 
    115          test = 1 
    116          IF (where(grille eq 'T'))[0] NE -1 THEN test = [test, gphit[*]] 
    117          IF (where(grille eq 'W'))[0] NE -1 AND (where(grille eq 'T'))[0] EQ -1 THEN $ 
    118           test = [test, gphit[*]] 
    119          IF (where(grille eq 'U'))[0] NE -1 AND n_elements(gphiu) NE 0 THEN test = [test, gphiu[*]] 
    120          IF (where(grille eq 'V'))[0] NE -1 AND n_elements(gphiv) NE 0 THEN test = [test, gphiv[*]] 
    121          IF (where(grille eq 'F'))[0] NE -1 AND n_elements(gphif) NE 0 THEN test = [test, gphif[*]] 
    122          test = test[1:n_elements(test)-1]  
    123          lat1=min([test], max = max) 
    124          lat2=max 
    125 ; prof1 et prof2 
    126          prof1=0 
    127          test = 1 
    128          IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN test = [test, gdept[jpk-1]] 
    129          IF (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN test = [test, gdepw[jpk-1]] 
    130          test = test[1:n_elements(test)-1]  
    131          prof2=max([test]) 
    132       end 
    133       1: BEGIN 
    134          xx1 = reform(x1) 
    135          taille = size(xx1) 
    136          CASE taille[1] OF 
    137             2:domdef, xx1[0], xx1[1]                                , FINDALWAYS = findalways,GRILLE=grille, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex 
    138             4:domdef, xx1[0], xx1[1], xx1[2], xx1[3]                , FINDALWAYS = findalways,GRILLE=grille, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex 
    139             6:domdef, xx1[0], xx1[1], xx1[2], xx1[3], xx1[4], xx1[5], FINDALWAYS = findalways,GRILLE=grille, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex 
    140             else: begin 
    141                ras = report('Mauvais nombre de parametre dans l''appel de DOMDEF') 
    142                return 
    143             end 
    144          ENDCASE 
    145          return 
    146       end 
    147       2: begin 
    148 ; lon1 et lon2 
    149          IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN BEGIN 
    150             ouca = where(gdept ge prof1 and gdept le prof2,nzt) 
    151             if ouca[0] NE -1 then zt=gdept[ouca] ELSE zt = -1 
    152          ENDIF 
    153          if (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 then BEGIN 
    154             ouca = where(gdepw ge prof1 and gdepw le prof2,nzw) 
    155             if ouca[0] NE -1 then zw=gdepw[ouca] ELSE zw = -1 
    156          ENDIF 
    157          test = 1 
    158          IF (where(grille eq 'T'))[0] NE -1 THEN test = [test, glamt[*]] 
    159          IF (where(grille eq 'W'))[0] NE -1 AND (where(grille eq 'T'))[0] EQ -1 THEN $ 
    160           test = [test, glamt[*]] 
    161          IF (where(grille eq 'U'))[0] NE -1 AND n_elements(glamu) NE 0 THEN test = [test, glamu[*]] 
    162          IF (where(grille eq 'V'))[0] NE -1 AND n_elements(glamv) NE 0 THEN test = [test, glamv[*]] 
    163          IF (where(grille eq 'F'))[0] NE -1 AND n_elements(glamf) NE 0 THEN test = [test, glamf[*]] 
    164          test = test[1:n_elements(test)-1]  
    165          lon1=min([test], max = max) 
    166          lon2=max 
    167 ; lat1 et lat2 
    168          test = 1 
    169          IF (where(grille eq 'T'))[0] NE -1 THEN test = [test, gphit[*]] 
    170          IF (where(grille eq 'W'))[0] NE -1 AND (where(grille eq 'T'))[0] EQ -1 THEN $ 
    171           test = [test, gphit[*]] 
    172          IF (where(grille eq 'U'))[0] NE -1 AND n_elements(gphiu) NE 0 THEN test = [test, gphiu[*]] 
    173          IF (where(grille eq 'V'))[0] NE -1 AND n_elements(gphiv) NE 0 THEN test = [test, gphiv[*]] 
    174          IF (where(grille eq 'F'))[0] NE -1 AND n_elements(gphif) NE 0 THEN test = [test, gphif[*]] 
    175          test = test[1:n_elements(test)-1]  
    176          lat1=min([test], max = max) 
    177          lat2=max 
    178 ; prof1 et prof2 
    179          if keyword_set(zindex) OR keyword_set(index) then BEGIN 
    180             z1 = x1 
    181             z2 = x2 
    182             IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN prof1 = gdept[z1] ELSE prof1=gdepw[z1] 
    183             if (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 then prof2=gdepw[z2] ELSE prof2=gdept[z2] 
    184          ENDIF ELSE BEGIN 
    185             prof1=x1 
    186             prof2=x2 
    187          ENDELSE  
    188       end 
    189       4: begin  
    190 ; lon1 et lon2 
    191          if keyword_set(xindex) OR keyword_set(index) then BEGIN 
    192             if keyword_set(yindex) OR keyword_set(index) then BEGIN 
    193                if n_elements(glamf) NE 0 then BEGIN 
    194                lon1 = min([glamt[x1:x2, y1:y2], glamf[x1:x2, y1:y2]], max = lon2) 
    195                ENDIF ELSE BEGIN  
    196                lon1 = min(glamt[x1:x2, y1:y2], max = lon2) 
    197                ENDELSE 
    198             ENDIF ELSE BEGIN  
    199                lon1 = glamt[x1, 0] 
    200                if n_elements(glamf) NE 0 then lon2 = glamf[x2, 0] $ 
    201                ELSE lon2 = glamt[x2, 0] 
     180; firstx[tuvf], lastx[tuvf], nx[tuvf] 
     181;------------------------------------------------------------------- 
     182; check if the grid is defined for U and V points. If not, take care 
     183; of the cases gridtype eq 'U' or 'V' 
     184;------------------------------------------------------------------- 
     185  errstatus = 0 
     186  IF (finite(glamu[0]*gphiu[0]) EQ 0 OR n_elements(glamu) EQ 0 OR n_elements(gphiu) EQ 0) AND (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN 
     187    firstxu = !values.f_nan 
     188    lastxu = !values.f_nan 
     189    nxu = !values.f_nan 
     190    okgrid = where(gridtype NE 'U', count) 
     191    IF count NE 0 THEN gridtype = gridtype[okgrid] $ 
     192    ELSE errstatus = report('U grid is undefined. Impossible to call domdef with vargid = ''U''') 
     193  ENDIF 
     194; 
     195  IF (finite(glamv[0]*gphiv[0]) EQ 0 OR n_elements(glamv) EQ 0 OR n_elements(gphiv) EQ 0) AND (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN 
     196    firstxv = !values.f_nan 
     197    lastxv = !values.f_nan 
     198    nxv = !values.f_nan 
     199    okgrid = where(gridtype NE 'V', count) 
     200    IF count NE 0 THEN gridtype = gridtype[okgrid] $ 
     201    ELSE errstatus = report('V grid is undefined. Impossible to call domdef with vargid = ''V''') 
     202  ENDIF 
     203  IF errstatus EQ -1 THEN return 
     204;------------------------------------------------------------------- 
     205;------------------------------------------------------------------- 
     206; horizontal domain defined with lon1, lon2, lat1 and lat2 
     207;------------------------------------------------------------------- 
     208;------------------------------------------------------------------- 
     209  IF N_PARAMS() EQ 0 $ 
     210    OR ( (N_PARAMS() EQ 4 OR N_PARAMS() EQ 6) $ 
     211         AND NOT keyword_set(xindex) AND NOT keyword_set(yindex) AND NOT keyword_set(index) ) THEN BEGIN 
     212    IF N_PARAMS() EQ 0 THEN BEGIN 
     213; find lon1 and lon2 the longitudinal boudaries of the full domain 
     214      IF (where(gridtype eq 'T'))[0] NE -1 THEN lon1t = min(glamt, max = lon2t) 
     215      IF (where(gridtype eq 'W'))[0] NE -1 AND (where(gridtype eq 'T'))[0] EQ -1 THEN lon1t = min(glamt, max = lon2t) 
     216      IF (where(gridtype eq 'U'))[0] NE -1 THEN lon1u = min(glamu, max = lon2u) 
     217      IF (where(gridtype eq 'V'))[0] NE -1 THEN lon1v = min(glamv, max = lon2v) 
     218      IF (where(gridtype eq 'F'))[0] NE -1 THEN lon1f = min(glamf, max = lon2f) 
     219      lon1 = min([lon1t, lon1u, lon1v, lon1f]) 
     220      lon2 = max([lon2t, lon2u, lon2v, lon2f]) 
     221; find lat1 and lat2 the latitudinal boudaries of the full domain 
     222      IF (where(gridtype eq 'T'))[0] NE -1 THEN lat1t = min(gphit, max = lat2t) 
     223      IF (where(gridtype eq 'W'))[0] NE -1 AND (where(gridtype eq 'T'))[0] EQ -1 THEN lat1t = min(gphit, max = lat2t) 
     224      IF (where(gridtype eq 'U'))[0] NE -1 THEN lat1u = min(gphiu, max = lat2u) 
     225      IF (where(gridtype eq 'V'))[0] NE -1 THEN lat1v = min(gphiv, max = lat2v) 
     226      IF (where(gridtype eq 'F'))[0] NE -1 THEN lat1f = min(gphif, max = lat2f) 
     227      lat1 = min([lat1t, lat1u, lat1v, lat1f]) 
     228      lat2 = max([lat2t, lat2u, lat2v, lat2f]) 
     229    ENDIF ELSE BEGIN  
     230      lon1 = min([x1, x2], max = lon2) 
     231      lat1 = min([y1, y2], max = lat2) 
     232    ENDELSE 
     233; find firstxt, firstxt, nxt and nyt according to lon1, lon2, lat1 and lat2 
     234    IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1 THEN BEGIN 
     235      dom = where( (glamt GE lon1) AND (glamt LE lon2) $ 
     236                   AND (gphit GE lat1) AND (gphit LE lat2) ) 
     237      IF (dom[0] EQ -1) THEN BEGIN 
     238        IF keyword_set(findalways) THEN BEGIN 
     239          domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 
     240          IF N_PARAMS() EQ 6 THEN $ 
     241            domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 
     242          RETURN 
     243        ENDIF ELSE BEGIN 
     244          ras = report('WARNING, The box does not contain any T points...') 
     245          firstxt = -1 & lastxt = -1 & nxt = 0 
     246          firstyt = -1 & lastyt = -1 & nyt = 0 
     247        ENDELSE 
     248      ENDIF ELSE BEGIN 
     249        jyt = dom / jpi 
     250        ixt = temporary(dom) MOD jpi 
     251        firstxt = min(temporary(ixt), max = lastxt) 
     252        firstyt = min(temporary(jyt), max = lastyt) 
     253        nxt = lastxt - firstxt + 1 
     254        nyt = lastyt - firstyt + 1 
     255      ENDELSE 
     256    ENDIF 
     257; find firstxu, firstxu, firstyu, firstyu, nxu and nyu 
     258; according to lon1, lon2, lat1 and lat2 
     259    IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN 
     260      IF keyword_set(memeindices) THEN BEGIN 
     261        firstxu = firstxt & lastxu = lastxt & nxu = nxt 
     262        firstyu = firstyt & lastyu = lastyt & nyu = nyt 
     263      ENDIF ELSE BEGIN 
     264        dom = where( (glamu GE lon1) AND (glamu LE lon2) $ 
     265                     AND (gphiu GE lat1) AND (gphiu LE lat2) ) 
     266        IF (dom[0] EQ -1) THEN BEGIN 
     267          IF keyword_set(findalways) THEN BEGIN 
     268            domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 
     269            IF N_PARAMS() EQ 6 THEN $ 
     270              domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 
     271            RETURN 
     272          ENDIF ELSE BEGIN 
     273            ras = report('WARNING, The box does not contain any U points...') 
     274            firstxu = -1 & lastxu = -1 & nxu = 0 
     275            firstyu = -1 & lastyu = -1 & nyu = 0 
     276          ENDELSE 
     277        ENDIF ELSE BEGIN 
     278          jyu = dom / jpi 
     279          ixu = temporary(dom) MOD jpi 
     280          firstxu = min(temporary(ixu), max = lastxu) 
     281          firstyu = min(temporary(jyu), max = lastyu) 
     282          nxu = lastxu - firstxu + 1 
     283          nyu = lastyu - firstyu + 1 
     284        ENDELSE 
     285      ENDELSE 
     286    ENDIF 
     287; find firstxv, firstxv, firstyv, firstyv, nxv and nyv 
     288;  according to lon1, lon2, lat1 and lat2 
     289    IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN 
     290      IF keyword_set(memeindices) THEN BEGIN 
     291        firstxv = firstxt & lastxv = lastxt & nxv = nxt 
     292        firstyv = firstyt & lastyv = lastyt & nyv = nyt 
     293      ENDIF ELSE BEGIN 
     294        dom = where( (glamv GE lon1) AND (glamv LE lon2) $ 
     295                     AND (gphiv GE lat1) AND (gphiv LE lat2) ) 
     296        IF (dom[0] EQ -1) THEN BEGIN 
     297          IF keyword_set(findalways) THEN BEGIN 
     298            domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 
     299            IF N_PARAMS() EQ 6 THEN $ 
     300              domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex  
     301            RETURN 
     302          ENDIF ELSE BEGIN 
     303            ras = report('WARNING, The box does not contain any V points...') 
     304            firstxv = -1 & lastxv = -1 & nxv = 0 
     305            firstyv = -1 & lastyv = -1 & nyv = 0 
     306          ENDELSE 
     307        ENDIF ELSE BEGIN 
     308          jyv = dom / jpi 
     309          ixv = temporary(dom) MOD jpi 
     310          firstxv = min(temporary(ixv), max = lastxv) 
     311          firstyv = min(temporary(jyv), max = lastyv) 
     312          nxv = lastxv - firstxv + 1 
     313          nyv = lastyv - firstyv + 1 
     314        ENDELSE 
     315      ENDELSE 
     316    ENDIF 
     317; find firstxf, firstxf, firstyf, firstyf, nxf and nyf 
     318; according to lon1, lon2, lat1 and lat2 
     319    IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN 
     320      IF keyword_set(memeindices) THEN BEGIN 
     321        firstxf = firstxt & lastxf = lastxt & nxf = nxt 
     322        firstyf = firstyt & lastyf = lastyt & nyf = nyt 
     323      ENDIF ELSE BEGIN 
     324        dom = where( (glamf GE lon1) AND (glamf LE lon2) $ 
     325                     AND (gphif GE lat1) AND (gphif LE lat2) ) 
     326        IF (dom[0] EQ -1) THEN BEGIN 
     327          IF keyword_set(findalways) THEN BEGIN 
     328            domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 
     329            IF N_PARAMS() EQ 6 THEN $ 
     330              domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 
     331            RETURN 
     332          ENDIF ELSE BEGIN 
     333            ras = report('WARNING, The box does not contain any F points...') 
     334            firstxf = -1 & lastxf = -1 & nxf = 0 
     335            firstyf = -1 & lastyf = -1 & nyf = 0 
     336          ENDELSE 
     337        ENDIF ELSE BEGIN 
     338          jyf = dom / jpi 
     339          ixf = temporary(dom) MOD jpi 
     340          firstxf = min(temporary(ixf), max = lastxf) 
     341          firstyf = min(temporary(jyf), max = lastyf) 
     342          nxf = lastxf - firstxf + 1 
     343          nyf = lastyf - firstyf + 1 
     344        ENDELSE 
     345      ENDELSE 
     346    ENDIF 
     347; 
     348  ENDIF ELSE BEGIN 
     349    CASE 1 OF 
     350;------------------------------------------------------------------- 
     351;------------------------------------------------------------------- 
     352; horizontal domain defined with the X and Y indexes 
     353;------------------------------------------------------------------- 
     354;------------------------------------------------------------------- 
     355      (keyword_set(xindex) AND keyword_set(yindex)) OR keyword_set(index):BEGIN 
     356        fstx = min([x1, x2], max = lstx) 
     357        fsty = min([y1, y2], max = lsty) 
     358        IF fstx LT 0 OR lstx GE jpi THEN BEGIN 
     359          ras = report('Bad definition of X1 or X2') 
     360          return 
     361        ENDIF 
     362        IF fsty LT 0 OR lsty GE jpj THEN BEGIN 
     363          ras = report('Bad definition of Y1 or Y2') 
     364          return 
     365        ENDIF 
     366        nx = lstx - fstx + 1 
     367        ny = lsty - fsty + 1 
     368; find lon1t, lon2t, lat1t, lat2t, firstxt, firstxt, nxt and nyt 
     369; according to x1, x2, y1, y2 
     370        IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype eq 'W'))[0] NE -1 THEN BEGIN 
     371          lon1t = min(glamt[fstx:lstx, fsty:lsty], max = lon2t) 
     372          lat1t = min(gphit[fstx:lstx, fsty:lsty], max = lat2t) 
     373          firstxt = fstx & lastxt = lstx 
     374          firstyt = fsty & lastyt = lsty 
     375          nxt = nx & nyt = ny 
     376        ENDIF 
     377; find lon1u, lon2u, lat1u, lat2u, firstxu, firstxu, nxu and nyu 
     378; according to x1, x2, y1, y2 
     379        IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN 
     380          lon1u = min(glamu[fstx:lstx, fsty:lsty], max = lon2u) 
     381          lat1u = min(gphiu[fstx:lstx, fsty:lsty], max = lat2u) 
     382          firstxu = fstx & lastxu = lstx 
     383          firstyu = fsty & lastyu = lsty 
     384          nxu = nx & nyu = ny 
     385        ENDIF 
     386; find lon1v, lon2v, lat1v, lat2v, firstxv, firstxv, nxv and nyv 
     387; according to x1, x2, y1, y2 
     388        IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN 
     389          lon1v = min(glamv[fstx:lstx, fsty:lsty], max = lon2v) 
     390          lat1v = min(gphiv[fstx:lstx, fsty:lsty], max = lat2v) 
     391          firstxv = fstx & lastxv = lstx 
     392          firstyv = fsty & lastyv = lsty 
     393          nxv = nx & nyv = ny 
     394        ENDIF          
     395; find lon1f, lon2f, lat1f, lat2f, firstxf, firstxf, nxf and nyf 
     396; according to x1, x2, y1, y2 
     397        IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN 
     398          lon1f = min(glamf[fstx:lstx, fsty:lsty], max = lon2f) 
     399          lat1f = min(gphif[fstx:lstx, fsty:lsty], max = lat2f) 
     400          firstxf = fstx & lastxf = lstx 
     401          firstyf = fsty & lastyf = lsty 
     402          nxf = nx & nyf = ny 
     403        ENDIF 
     404        lon1 = min([lon1t, lon1u, lon1v, lon1f]) 
     405        lon2 = max([lon2t, lon2u, lon2v, lon2f]) 
     406        lat1 = min([lat1t, lat1u, lat1v, lat1f]) 
     407        lat2 = max([lat2t, lat2u, lat2v, lat2f]) 
     408      END 
     409;------------------------------------------------------------------- 
     410;------------------------------------------------------------------- 
     411; horizontal domain defined with the X index and lat1/lat2 
     412;------------------------------------------------------------------- 
     413;------------------------------------------------------------------- 
     414      keyword_set(xindex):BEGIN 
     415        fstx = min([x1, x2], max = lstx) 
     416        IF fstx LT 0 OR lstx GE jpi THEN BEGIN 
     417          ras = report('Bad definition of X1 or X2') 
     418          return 
     419        ENDIF 
     420        nx = lstx - fstx + 1 
     421        lat1 = min([y1, y2], max = lat2) 
     422; find lon1t, lon2t, firstxt, firstxt, firstyt, firstyt, nxt 
     423; and nyt according to x1, x2, lat1 and lat2 
     424        IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1 THEN BEGIN 
     425          firstxt = fstx & lastxt = lstx & nxt = nx 
     426          dom = where( (gphit[fstx:lstx, *] GE lat1) AND (gphit[fstx:lstx, *] LE lat2) ) 
     427          IF (dom[0] EQ -1) THEN BEGIN 
     428            IF keyword_set(findalways) THEN BEGIN 
     429              domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 
     430              IF N_PARAMS() EQ 6 THEN $ 
     431                domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 
     432              RETURN 
     433            ENDIF ELSE BEGIN 
     434              ras = report('WARNING, The box does not contain any T points...') 
     435              firstyt = -1 & lastyt = -1 & nyt = 0 
    202436            ENDELSE 
    203          ENDIF ELSE BEGIN 
    204             lon1=x1 
    205             lon2=x2 
    206          ENDELSE  
    207 ; lat1 et lat2 
    208          if keyword_set(yindex) OR keyword_set(index) then BEGIN 
    209             if keyword_set(xindex) OR keyword_set(index) then begin 
    210                lat1 = min([gphit[x1, y1], gphit[x2, y1]]) 
    211                if n_elements(gphif) NE 0 then lat2 = max([gphif[x1, y2], gphif[x2, y2]]) $ 
    212                ELSE lat2 = max([gphit[x1, y2], gphit[x2, y2]]) 
    213             ENDIF ELSE BEGIN  
    214                lat1 = gphit[0, y1] 
    215                if n_elements(gphif) NE 0 then lat2 = gphif[0, y2] $ 
    216                ELSE lat2 = gphit[0, y2] 
     437          ENDIF ELSE BEGIN 
     438            jyt = temporary(dom) / nx 
     439            firstyt = min(temporary(jyt), max = lastyt) 
     440            nyt = lastyt - firstyt + 1 
     441          ENDELSE 
     442          IF nyt NE 0 THEN lon1t = min(glamt[firstxt:lastxt, firstyt:lastyt], max = lon2t) 
     443        ENDIF 
     444; find lon1u, lon2u, firstxu, firstxu, firstyu, firstyu, nxu 
     445; and nyu according to x1, x2, lat1 and lat2 
     446        IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN 
     447          firstxu = fstx & lastxu = lstx & nxu = nx 
     448          IF keyword_set(memeindices) THEN BEGIN 
     449            firstyu = firstyt & lastyu = lastyt & nyu = nyt 
     450          ENDIF ELSE BEGIN 
     451            dom = where( (gphiu[fstx:lstx, *] GE lat1) AND (gphiu[fstx:lstx, *] LE lat2) ) 
     452            IF (dom[0] EQ -1) THEN BEGIN 
     453              IF keyword_set(findalways) THEN BEGIN 
     454                domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 
     455                IF N_PARAMS() EQ 6 THEN $ 
     456                  domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 
     457                RETURN 
     458              ENDIF ELSE BEGIN 
     459                ras = report('WARNING, The box does not contain any U points...') 
     460                firstyu = -1 & lastyu = -1 & nyu = 0 
     461              ENDELSE 
     462            ENDIF ELSE BEGIN 
     463              jyu = temporary(dom) / nx 
     464              firstyu = min(temporary(jyu), max = lastyu) 
     465              nyu = lastyu - firstyu + 1 
    217466            ENDELSE 
    218          ENDIF ELSE BEGIN 
    219             lat1=y1 
    220             lat2=y2 
    221          ENDELSE  
    222 ; prof1 et prof2 
    223          prof1=0 
    224          test = 1 
    225          IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN test = [test, gdept[jpk-1]] 
    226          IF (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN test = [test, gdepw[jpk-1]] 
    227          test = test[1:n_elements(test)-1]  
    228          prof2=max([test]) 
    229       end 
    230       6: begin  
    231 ; lon1 et lon2 
    232          if keyword_set(xindex) OR keyword_set(index) then BEGIN 
    233             if keyword_set(yindex) OR keyword_set(index) then begin 
    234                lon1 = min([glamt[x1, y1], glamt[x1, y2]]) 
    235                if n_elements(glamf) NE 0 then lon2 = max([glamf[x2, y1], glamf[x2, y2]]) $ 
    236                ELSE lon2 = max([glamt[x2, y1], glamt[x2, y2]]) 
    237             ENDIF ELSE BEGIN  
    238                lon1 = glamt[x1, 0] 
    239                if n_elements(glamf) NE 0 then lon2 = glamf[x2, 0] $ 
    240                ELSE lon2 = glamt[x2, 0] 
     467          ENDELSE 
     468          IF nyu NE 0 THEN lon1u = min(glamu[firstxu:lastxu, firstyu:lastyu], max = lon2u) 
     469        ENDIF 
     470; find lon1v, lon2v, firstxv, firstxv, firstyv, firstyv, nxv 
     471; and nyv according to x1, x2, lat1 and lat2 
     472        IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN 
     473          firstxv = fstx & lastxv = lstx & nxv = nx 
     474          IF keyword_set(memeindices) THEN BEGIN 
     475            firstyv = firstyt & lastyv = lastyt & nyv = nyt 
     476          ENDIF ELSE BEGIN 
     477            dom = where( (gphiv[fstx:lstx, *] GE lat1) AND (gphiv[fstx:lstx, *] LE lat2) ) 
     478            IF (dom[0] EQ -1) THEN BEGIN 
     479              IF keyword_set(findalways) THEN BEGIN 
     480                domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 
     481                IF N_PARAMS() EQ 6 THEN $ 
     482                  domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 
     483                RETURN 
     484              ENDIF ELSE BEGIN 
     485                ras = report('WARNING, The box does not contain any V points...') 
     486                firstyv = -1 & lastyv = -1 & nyv = 0 
     487              ENDELSE 
     488            ENDIF ELSE BEGIN 
     489              jyv = temporary(dom) / nx 
     490              firstyv = min(temporary(jyv), max = lastyv) 
     491              nyv = lastyv - firstyv + 1 
    241492            ENDELSE 
    242          ENDIF ELSE BEGIN 
    243             lon1=x1 
    244             lon2=x2 
    245          ENDELSE  
    246 ; lat1 et lat2 
    247          if keyword_set(yindex) OR keyword_set(index) then BEGIN 
    248             if keyword_set(xindex) OR keyword_set(index) then begin 
    249                lat1 = min([gphit[x1, y1], gphit[x2, y1]]) 
    250                if n_elements(gphif) NE 0 then lat2 = max([gphif[x1, y2], gphif[x2, y2]]) $ 
    251                ELSE lat2 = max([gphit[x1, y2], gphit[x2, y2]]) 
    252             ENDIF ELSE BEGIN  
    253                lat1 = gphit[0, y1] 
    254                if n_elements(gphif) NE 0 then lat2 = gphif[0, y2] $ 
    255                ELSE lat2 = gphit[0, y2] 
     493          ENDELSE 
     494          IF nyv NE 0 THEN lon1v = min(glamv[firstxv:lastxv, firstyv:lastyv], max = lon2v) 
     495        ENDIF 
     496; find lon1f, lon2f, firstxf, firstxf, firstyf, firstyf, nxf 
     497; and nyf according to x1, x2, lat1 and lat2 
     498        IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN 
     499          firstxf = fstx & lastxf = lstx & nxf = nx 
     500          IF keyword_set(memeindices) THEN BEGIN 
     501            firstyf = firstyt & lastyf = lastyt & nyf = nyt 
     502          ENDIF ELSE BEGIN 
     503            dom = where( (gphif[fstx:lstx, *] GE lat1) AND (gphif[fstx:lstx, *] LE lat2) ) 
     504            IF (dom[0] EQ -1) THEN BEGIN 
     505              IF keyword_set(findalways) THEN BEGIN 
     506                domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 
     507                IF N_PARAMS() EQ 6 THEN $ 
     508                  domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 
     509                RETURN 
     510              ENDIF ELSE BEGIN 
     511                ras = report('WARNING, The box does not contain any F points...') 
     512                firstyf = -1 & lastyf = -1 & nyf = 0 
     513              ENDELSE 
     514            ENDIF ELSE BEGIN 
     515              jyf = temporary(dom) / nx 
     516              firstyf = min(temporary(jyf), max = lastyf) 
     517              nyf = lastyf - firstyf + 1 
    256518            ENDELSE 
    257          ENDIF ELSE BEGIN 
    258             lat1=y1 
    259             lat2=y2 
    260          ENDELSE  
    261 ; prof1 et prof2 
    262          if keyword_set(zindex) OR keyword_set(index) then begin 
    263             IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN prof1=gdept[z1] ELSE prof1=gdepw[z1] 
    264             if (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 then prof2=gdepw[z2] ELSE prof2=gdept[z2] 
    265          ENDIF ELSE BEGIN 
    266             prof1=z1 
    267             prof2=z2 
    268          ENDELSE  
    269       end 
    270       else: begin 
    271          ras = report('Mauvais nombre de parametre dans l''appel de DOMDEF') 
    272          return 
    273       end 
    274    endcase   
    275 ;----------------------------------------------------------------------        
    276 ; determination de nzt et nwt 
    277 ;----------------------------------------------------------------------        
    278    IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN BEGIN 
    279       if keyword_set(zindex) OR keyword_set(index) then nzt = z2-z1+1 $ 
    280       ELSE BEGIN  
    281          ouca = where(gdept ge prof1 and gdept le prof2,nzt) 
    282          if ouca[0] NE -1 then zt=gdept[ouca] ELSE zt = -1 
    283       ENDELSE 
    284    ENDIF 
    285    if (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 then BEGIN 
    286       if keyword_set(zindex) OR keyword_set(index) then nzw = z2-z1+1 $ 
    287       ELSE BEGIN  
    288          ouca = where(gdepw ge prof1 and gdepw le prof2,nzw) 
    289          if ouca[0] NE -1 then zw=gdepw[ouca] ELSE zw = -1 
    290       ENDELSE 
    291    ENDIF 
    292 ;----------------------------------------------------------------------        
    293 ; determination de premierzt, dernierzt, premierzw, dernierzw 
    294 ;----------------------------------------------------------------------        
    295    IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN BEGIN 
    296       if keyword_set(zindex) OR keyword_set(index) then BEGIN  
    297          premierzt = z1 
    298          dernierzt = z2 
     519          ENDELSE 
     520          IF nyf NE 0 THEN lon1f = min(glamf[firstxf:lastxf, firstyf:lastyf], max = lon2f) 
     521        ENDIF 
     522        lon1 = min([lon1t, lon1u, lon1v, lon1f]) 
     523        lon2 = max([lon2t, lon2u, lon2v, lon2f]) 
     524      END 
     525;------------------------------------------------------------------- 
     526;------------------------------------------------------------------- 
     527; horizontal domain defined with the Y index and lon1/lon2 
     528;------------------------------------------------------------------- 
     529;------------------------------------------------------------------- 
     530      keyword_set(yindex):BEGIN 
     531        fsty = min([y1, y2], max = lsty) 
     532        IF fsty LT 0 OR lsty GE jpj THEN BEGIN 
     533          ras = report('Bad definition of Y1 or Y2') 
     534          return 
     535        ENDIF 
     536        ny = lsty - fsty + 1 
     537        lon1 = min([x1, x2], max = lon2) 
     538; find lat1t, lat2t, firstxt, firstxt, firstyt, firstyt, nxt 
     539; and nyt according to x1, x2, lon1 and lon2 
     540        IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1 THEN BEGIN 
     541          firstyt = fsty & lastyt = lsty & nyt = ny 
     542          dom = where( (glamt[*, fsty:lsty] GE lon1) AND (glamt[*, fsty:lsty] LE lon2) ) 
     543          IF (dom[0] EQ -1) THEN BEGIN 
     544            IF keyword_set(findalways) THEN BEGIN 
     545              domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 
     546              IF N_PARAMS() EQ 6 THEN $ 
     547                domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex  
     548              RETURN 
     549            ENDIF ELSE BEGIN 
     550              ras = report('WARNING, The box does not contain any T points...') 
     551              firstxt = -1 & lastxt = -1 & nxt = 0 
     552            ENDELSE 
     553          ENDIF ELSE BEGIN 
     554            jxt = temporary(dom) MOD jpi 
     555            firstxt = min(temporary(jxt), max = lastxt) 
     556            nxt = lastxt - firstxt + 1 
     557          ENDELSE 
     558          IF nxt NE 0 THEN lat1t = min(gphit[firstxt:lastxt, firstyt:lastyt], max = lat2t) 
     559        ENDIF 
     560; find lat1u, lat2u, firstxu, firstxu, firstyu, firstyu, nxu 
     561; and nyu according to x1, x2, lon1 and lon2 
     562        IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN 
     563          firstyu = fsty & lastyu = lsty & nyu = ny 
     564          IF keyword_set(memeindices) THEN BEGIN 
     565            firstxu = firstyt & lastxu = lastyt & nxu = nxt 
     566          ENDIF ELSE BEGIN 
     567            dom = where( (glamu[*, fsty:lsty] GE lon1) AND (glamu[*, fsty:lsty] LE lon2) ) 
     568            IF (dom[0] EQ -1) THEN BEGIN 
     569              IF keyword_set(findalways) THEN BEGIN 
     570                domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 
     571                IF N_PARAMS() EQ 6 THEN $ 
     572                  domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex  
     573                RETURN 
     574              ENDIF ELSE BEGIN 
     575                ras = report('WARNING, The box does not contain any U points...') 
     576                firstxu = -1 & lastxu = -1 & nxu = 0 
     577              ENDELSE 
     578            ENDIF ELSE BEGIN 
     579              jxu = temporary(dom) MOD jpi 
     580              firstxu = min(temporary(jxu), max = lastxu) 
     581              nxu = lastxu - firstxu + 1 
     582            ENDELSE 
     583          ENDELSE 
     584          IF nxu NE 0 THEN lat1u = min(gphiu[firstxu:lastxu, firstyu:lastyu], max = lat2u) 
     585        ENDIF 
     586; find lat1v, lat2v, firstxv, firstxv, firstyv, firstyv, nxv 
     587; and nyv according to x1, x2, lon1 and lon2 
     588        IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN 
     589          firstyv = fsty & lastyv = lsty & nyv = ny 
     590          IF keyword_set(memeindices) THEN BEGIN 
     591            firstxv = firstyt & lastxv = lastyt & nxv = nxt 
     592          ENDIF ELSE BEGIN 
     593            dom = where( (glamv[*, fsty:lsty] GE lon1) AND (glamv[*, fsty:lsty] LE lon2) ) 
     594            IF (dom[0] EQ -1) THEN BEGIN 
     595              IF keyword_set(findalways) THEN BEGIN 
     596                domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 
     597                IF N_PARAMS() EQ 6 THEN $ 
     598                  domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex  
     599                RETURN 
     600              ENDIF ELSE BEGIN 
     601                ras = report('WARNING, The box does not contain any V points...') 
     602                firstxv = -1 & lastxv = -1 & nxv = 0 
     603              ENDELSE 
     604            ENDIF ELSE BEGIN 
     605              jxv = temporary(dom) MOD jpi 
     606              firstxv = min(temporary(jxv), max = lastxv) 
     607              nxv = lastxv - firstxv + 1 
     608            ENDELSE 
     609          ENDELSE 
     610          IF nxv NE 0 THEN lat1v = min(gphiv[firstxv:lastxv, firstyv:lastyv], max = lat2v) 
     611        ENDIF 
     612; find lat1f, lat2f, firstxf, firstxf, firstyf, firstyf, nxf 
     613; and nyf according to x1, x2, lon1 and lon2 
     614        IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN 
     615          firstyf = fsty & lastyf = lsty & nyf = ny 
     616          IF keyword_set(memeindices) THEN BEGIN 
     617            firstxf = firstyt & lastxf = lastyt & nxf = nxt 
     618          ENDIF ELSE BEGIN 
     619            dom = where( (glamf[*, fsty:lsty] GE lon1) AND (glamf[*, fsty:lsty] LE lon2) ) 
     620            IF (dom[0] EQ -1) THEN BEGIN 
     621              IF keyword_set(findalways) THEN BEGIN 
     622                domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex 
     623                IF N_PARAMS() EQ 6 THEN $ 
     624                  domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex  
     625                RETURN 
     626              ENDIF ELSE BEGIN 
     627                ras = report('WARNING, The box does not contain any F points...') 
     628                firstxf = -1 & lastyf = -1 & nxf = 0 
     629              ENDELSE 
     630            ENDIF ELSE BEGIN 
     631              jxf = temporary(dom) MOD jpi 
     632              firstxf = min(temporary(jxf), max = lastxf) 
     633              nxf = lastxf - firstxf + 1 
     634            ENDELSE 
     635          ENDELSE 
     636          IF nxf NE 0 THEN lat1f = min(gphif[firstxf:lastxf, firstyf:lastyf], max = lat2f) 
     637        ENDIF 
     638        lat1 = min([lat1t, lat1u, lat1v, lat1f]) 
     639        lat2 = max([lat2t, lat2u, lat2v, lat2f]) 
     640      END 
     641    ENDCASE 
     642  ENDELSE  
     643;------------------------------------------------------------------- 
     644; The extracted domain is it regular or not?  
     645;------------------------------------------------------------------- 
     646  CASE 1 OF 
     647    ((where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype eq 'W'))[0] NE -1) AND nxt NE 0 AND nyt NE 0:BEGIN 
     648; to get faster, we first test the most basic cases befor testing the 
     649; full array. 
     650      CASE 0 OF 
     651        array_equal(glamt[firstxt:lastxt, firstyt], glamt[firstxt:lastxt, lastyt]):key_irregular = 1b 
     652        array_equal(gphit[firstxt, firstyt:lastyt], gphit[lastxt, firstyt:lastyt]):key_irregular = 1b 
     653        array_equal(glamt[firstxt:lastxt, firstyt:lastyt] $ 
     654                    , glamt[firstxt:lastxt, firstyt]#replicate(1, nyt)):key_irregular = 1b 
     655        array_equal(gphit[firstxt:lastxt, firstyt:lastyt] $ 
     656                    , replicate(1, nxt)#(gphit[firstxt, firstyt:lastyt])[*]):key_irregular = 1b 
     657        ELSE:key_irregular = 0b 
     658      ENDCASE 
     659    END 
     660    (where(gridtype eq 'U'))[0] NE -1 AND nxu NE 0 AND nyu NE 0:BEGIN 
     661      CASE 0 OF 
     662        array_equal(glamu[firstxu:lastxu, firstyu], glamu[firstxu:lastxu, lastyu]):key_irregular = 1b 
     663        array_equal(gphiu[firstxu, firstyu:lastyu], gphiu[lastxu, firstyu:lastyu]):key_irregular = 1b 
     664        array_equal(glamu[firstxu:lastxu, firstyu:lastyu] $ 
     665                    , glamu[firstxu:lastxu, firstyu]#replicate(1, nyu)):key_irregular = 1b 
     666        array_equal(gphiu[firstxu:lastxu, firstyu:lastyu] $ 
     667                    , replicate(1, nxu)#(gphiu[firstxu, firstyu:lastyu])[*]):key_irregular = 1b 
     668        ELSE:key_irregular = 0b 
     669      ENDCASE 
     670    END 
     671    (where(gridtype eq 'V'))[0] NE -1 AND nxv NE 0 AND nyv NE 0:BEGIN 
     672      CASE 0 OF 
     673        array_equal(glamv[firstxv:lastxv, firstyv], glamv[firstxv:lastxv, lastyv]):key_irregular = 1b 
     674        array_equal(gphiv[firstxv, firstyv:lastyv], gphiv[lastxv, firstyv:lastyv]):key_irregular = 1b 
     675        array_equal(glamv[firstxv:lastxv, firstyv:lastyv] $ 
     676                    , glamv[firstxv:lastxv, firstyv]#replicate(1, nyv)):key_irregular = 1b 
     677        array_equal(gphiv[firstxv:lastxv, firstyv:lastyv] $ 
     678                    , replicate(1, nxv)#(gphiv[firstxv, firstyv:lastyv])[*]):key_irregular = 1b 
     679        ELSE:key_irregular = 0b 
     680      ENDCASE 
     681    END 
     682    (where(gridtype eq 'F'))[0] NE -1 AND nxf NE 0 AND nyf NE 0:BEGIN 
     683      CASE 0 OF 
     684        array_equal(glamf[firstxf:lastxf, firstyf], glamf[firstxf:lastxf, lastyf]):key_irregular = 1b 
     685        array_equal(gphif[firstxf, firstyf:lastyf], gphif[lastxf, firstyf:lastyf]):key_irregular = 1b 
     686        array_equal(glamf[firstxf:lastxf, firstyf:lastyf] $ 
     687                    , glamf[firstxf:lastxf, firstyf]#replicate(1, nyf)):key_irregular = 1b 
     688        array_equal(gphif[firstxf:lastxf, firstyf:lastyf] $ 
     689                    , replicate(1, nxf)#(gphif[firstxf, firstyf:lastyf])[*]):key_irregular = 1b 
     690        ELSE:key_irregular = 0b 
     691      ENDCASE 
     692    END 
     693    ELSE: 
     694  ENDCASE 
     695; 
     696;------------------------------------------------------------------- 
     697;------------------------------------------------------------------- 
     698; define all vertical parameters ... 
     699; vert1, vert2 
     700; firstz[tw], lastz[tw], nz[tw] 
     701;------------------------------------------------------------------- 
     702;------------------------------------------------------------------- 
     703; 
     704  vertical: 
     705; 
     706;------------------------------------------------------------------- 
     707; vertical domain defined with vert1, vert2 
     708;------------------------------------------------------------------- 
     709  IF NOT (keyword_set(zindex) OR keyword_set(index)) THEN BEGIN 
     710; define vert1 et vert2 
     711    CASE N_PARAMS() OF 
     712      2:vert1 = min([x1, x2], max = vert2) 
     713      6:vert1 = min([z1, z2], max = vert2) 
     714      ELSE:BEGIN 
     715        IF (inter(byte(gridtype), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN $ 
     716          vert1t = min(gdept, max = vert2t) 
     717        IF (where(gridtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN $ 
     718          vert1w = min(gdepw, max = vert2w) 
     719        vert1 = min([vert1t, vert1w]) 
     720        vert2 = max([vert2t, vert2w]) 
     721      END 
     722    ENDCASE 
     723; define firstzt, firstzt, nzt 
     724    IF (inter(byte(gridtype), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN BEGIN 
     725      domz = where(gdept ge vert1 and gdept le vert2, nzt) 
     726      IF nzt NE 0 THEN BEGIN 
     727        firstzt = domz[0] 
     728        lastzt = domz[nzt-1] 
    299729      ENDIF ELSE BEGIN  
    300          if zt[0] NE -1 then begin 
    301             premierzt = (where(gdept eq zt[0]))[0] 
    302             dernierzt = (where(gdept eq zt[nzt-1]))[0] 
    303          ENDIF ELSE BEGIN  
    304             premierzt = -1 
    305             dernierzt = -1 
    306          ENDELSE 
    307       ENDELSE 
    308    endif 
    309    if (where(grille eq 'W'))[0] NE -1  AND n_elements(gdepw) NE 0 then begin 
    310       if keyword_set(zindex) OR keyword_set(index) then BEGIN  
    311          premierzw = z1 
    312          dernierzw = z2 
     730        ras = report('WARNING, The box does not contain any T level...') 
     731        firstzt = -1 
     732        lastzt = -1 
     733      ENDELSE  
     734    ENDIF 
     735; define firstzw, firstzw, nzw 
     736    IF (where(gridtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN BEGIN 
     737      IF keyword_set(memeindices) THEN BEGIN 
     738        firstzw = firstzt 
     739        lastzw = lastzt 
     740        nzw = nzt 
    313741      ENDIF ELSE BEGIN  
    314          if zw[0] NE -1 then begin 
    315             premierzw = (where(gdepw eq zw[0]))[0] 
    316             dernierzw = (where(gdepw eq zw[nzw-1]))[0] 
    317          ENDIF ELSE BEGIN  
    318             premierzw = -1 
    319             dernierzw = -1 
    320          ENDELSE 
    321       ENDELSE 
    322    endif 
    323 ;------------------------------------------------------------------- 
    324 ;------------------------------------------------------------------- 
    325 ; calcul pour les pts t et w 
    326 ;------------------------------------------------------------------- 
    327 ;------------------------------------------------------------------- 
    328    tempdeux = systime(1)        ; pour key_performance =2 
    329    if (where(grille eq 'T'))[0] NE -1 or (where(grille EQ 'W'))[0] NE -1 then begin 
    330       if (keyword_set(xindex) AND keyword_set(yindex)) OR keyword_set(index) then begin 
    331          premierxt = x1 
    332          premieryt = y1 
    333          dernierxt = x2 
    334          dernieryt = y2 
    335       ENDIF ELSE BEGIN 
    336          dom =   where( (glamt ge lon1) $ 
    337                         and (glamt le lon2) $ 
    338                         and (gphit ge lat1) $ 
    339                         and (gphit le lat2) ) 
    340          if (dom[0] eq -1) then BEGIN 
    341             if keyword_set(findalways) then BEGIN 
    342                domdef, prof1,prof2, GRILLE=grille, MEMEINDICES = memeindices, _extra = ex 
    343                return 
    344             ENDIF ELSE BEGIN 
    345                ras = report('ATTENTION, la boite ne contient pas de points t') 
    346                goto, ptu 
    347             ENDELSE 
    348          endif  
    349          IF testvar(var = key_performance) EQ 2 THEN $ 
    350           print, 'temps domdef: trouver le sous domaine', systime(1)-tempdeux 
    351 ;------------------------------------------------------------------- 
    352 ; on recupere les abscisses et les ordonees des points selectionnes 
    353 ;------------------------------------------------------------------- 
    354          tempdeux = systime(1)  ; pour key_performance =2 
    355          ordonnee = dom/jpi 
    356          abscisse = dom-ordonnee*jpi 
    357 ;------------------------------------------------------------------- 
    358 ; on les ordonnes et on elimine les elements en double 
    359 ;------------------------------------------------------------------- 
    360          ordonnee = ordonnee 
    361          abscisse = abscisse 
    362          ordonnee = ordonnee(uniq(ordonnee, sort(ordonnee))) 
    363          abscisse = abscisse(uniq(abscisse, sort(abscisse))) 
    364 ;------------------------------------------------------------------- 
    365 ; indices des bornes du domaine 
    366 ;------------------------------------------------------------------- 
    367          premierxt = abscisse[0] 
    368          premieryt = ordonnee[0] 
    369          dernierxt = abscisse[n_elements(abscisse) -1] 
    370          dernieryt = ordonnee[n_elements(ordonnee) -1] 
     742        domz = where(gdepw ge vert1 and gdepw le vert2, nzw) 
     743        IF nzw NE 0 THEN BEGIN 
     744          firstzw = domz[0] 
     745          lastzw = domz[nzw-1] 
     746        ENDIF ELSE BEGIN  
     747          ras = report('WARNING, The box does not contain any W level...') 
     748          firstzw = -1 
     749          lastzw = -1 
     750        ENDELSE  
    371751      ENDELSE  
    372 ;------------------------------------------------------------------- 
    373 ; nombre de points selectionnes 
    374 ;------------------------------------------------------------------- 
    375       nxt = dernierxt-premierxt+1 
    376       nyt = dernieryt-premieryt+1 
    377       IF testvar(var = key_performance) EQ 2 THEN $ 
    378        print, 'temps domdef: calculs de premier, dernier,...', systime(1)-tempdeux 
    379 ;------------------------------------------------------------------- 
    380 ; la grille est reguliere ou non ?? 
    381 ;------------------------------------------------------------------- 
    382       tempdeux = systime(1)     ; pour key_performance =2 
    383       if (total(glamt[premierxt:dernierxt,premieryt:dernieryt] $ 
    384                 NE glamt[premierxt:dernierxt,0]#replicate(1, nyt)) eq 0 $ 
    385           AND total(gphit[premierxt:dernierxt,premieryt:dernieryt] $ 
    386                     NE replicate(1, nxt)#(gphit[0,premieryt:dernieryt])[*]) eq 0 ) $ 
    387        then key_irregular = 0 ELSE key_irregular = 1 
    388       IF testvar(var = key_performance) EQ 2 THEN $ 
    389        print, 'temps domdef: la grille est reguliere ou non ', systime(1)-tempdeux 
    390 ;------------------------------------------------------------------- 
    391       if keyword_set(memeindices) then begin 
    392          premierxu = premierxt & premierxv = premierxt & premierxf = premierxt  
    393          dernierxu = dernierxt & dernierxv = dernierxt & dernierxf = dernierxt 
    394          premieryu = premieryt & premieryv = premieryt & premieryf = premieryt  
    395          dernieryu = dernieryt & dernieryv = dernieryt & dernieryf = dernieryt 
    396          premierzw = premierzt  
    397          dernierzw = dernierzt  
    398          nxu = nxt & nxv = nxt & nxf = nxt 
    399          nyu = nyt & nyv = nyt & nyf = nyt 
    400          nzw = nzt 
    401          lon1 = glamt[premierxt, 0] 
    402          lon2 = glamf[dernierxf, 0] 
    403          lat1 = gphit[0, premieryt] 
    404          lat2 = gphif[0, dernieryf] 
    405          return 
    406       endif 
    407 ptu: 
    408    endif 
    409 ;------------------------------------------------------------------- 
    410 ; calcul pour les pts u 
    411 ;------------------------------------------------------------------- 
    412    if (where(grille eq 'U'))[0] NE -1 then begin 
    413       if (keyword_set(xindex) AND keyword_set(xindex)) OR keyword_set(index) then begin 
    414          premierxu = x1 
    415          premieryu = y1 
    416          dernierxu = x2 
    417          dernieryu = y2 
    418       ENDIF ELSE BEGIN 
    419          dom =   where( (glamu ge lon1) $ 
    420                         and (glamu le lon2) $ 
    421                         and (gphiu ge lat1) $ 
    422                         and (gphiu le lat2) ) 
    423          if (dom[0] eq -1) then begin 
    424             if keyword_set(findalways) then begin 
    425                domdef, prof1,prof2, GRILLE=grille, MEMEINDICES = memeindices, _extra = ex 
    426                return 
    427             ENDIF ELSE BEGIN 
    428                ras = report('ATTENTION, la boite ne contient pas de points u') 
    429                goto, ptv 
    430             ENDELSE 
    431          ENDIF 
    432          ordonnee = dom/jpi 
    433          abscisse = dom-ordonnee*jpi 
    434          ordonnee = ordonnee(uniq(ordonnee, sort(ordonnee))) 
    435          abscisse = abscisse(uniq(abscisse, sort(abscisse))) 
    436          premierxu = abscisse[0] 
    437          premieryu = ordonnee[0] 
    438          dernierxu = abscisse[n_elements(abscisse) -1] 
    439          dernieryu = ordonnee[n_elements(ordonnee) -1] 
    440       ENDELSE  
    441       nxu = dernierxu-premierxu+1 
    442       nyu = dernieryu-premieryu+1 
    443       if (total(glamu[premierxu:dernierxu,premieryu:dernieryu] $ 
    444                 NE glamu[premierxu:dernierxu,0]#replicate(1, nyu)) eq 0 $ 
    445           AND total(gphiu[premierxu:dernierxu,premieryu:dernieryu] $ 
    446                     NE replicate(1, nxu)#(gphiu[0,premieryu:dernieryu])[*]) eq 0 ) $ 
    447        then key_irregular = 0 ELSE key_irregular = 1 
    448 ptv: 
    449    endif 
    450 ;------------------------------------------------------------------- 
    451 ; calcul pour les pts v 
    452 ;------------------------------------------------------------------- 
    453    if (where(grille eq 'V'))[0] NE -1 then begin 
    454       if (keyword_set(xindex) AND keyword_set(xindex)) OR keyword_set(index) then begin 
    455          premierxv = x1 
    456          premieryv = y1 
    457          dernierxv = x2 
    458          dernieryv = y2 
    459       ENDIF ELSE BEGIN 
    460          dom =   where( (glamv ge lon1) $ 
    461                         and (glamv le lon2) $ 
    462                         and (gphiv ge lat1) $ 
    463                         and (gphiv le lat2) ) 
    464          if (dom[0] eq -1) then begin 
    465             if keyword_set(findalways) then begin 
    466                domdef, prof1,prof2, GRILLE=grille, MEMEINDICES = memeindices, _extra = ex 
    467                return 
    468             ENDIF ELSE BEGIN 
    469                ras = report('ATTENTION, la boite ne contient pas de points v') 
    470                goto, ptf 
    471             ENDELSE 
    472          ENDIF 
    473          ordonnee = dom/jpi 
    474          abscisse = dom-ordonnee*jpi 
    475          ordonnee = ordonnee(uniq(ordonnee, sort(ordonnee))) 
    476          abscisse = abscisse(uniq(abscisse, sort(abscisse))) 
    477          premierxv = abscisse[0] 
    478          premieryv = ordonnee[0] 
    479          dernierxv = abscisse[n_elements(abscisse) -1] 
    480          dernieryv = ordonnee[n_elements(ordonnee) -1] 
    481       ENDELSE  
    482       nxv = dernierxv-premierxv+1 
    483       nyv = dernieryv-premieryv+1 
    484       if (total(glamv[premierxv:dernierxv,premieryv:dernieryv] $ 
    485                 NE glamv[premierxv:dernierxv,0]#replicate(1, nyv)) eq 0 $ 
    486           AND total(gphiv[premierxv:dernierxv,premieryv:dernieryv] $ 
    487                     NE replicate(1, nxv)#(gphiv[0,premieryv:dernieryv])[*]) eq 0 ) $ 
    488        then key_irregular = 0 ELSE key_irregular = 1 
    489 ptf: 
    490    endif 
    491 ;------------------------------------------------------------------- 
    492 ; calcul pour les pts f 
    493 ;------------------------------------------------------------------- 
    494    if (where(grille eq 'F'))[0] NE -1 then begin 
    495       if (keyword_set(xindex) AND keyword_set(xindex)) OR keyword_set(index) then begin 
    496          premierxf = x1 
    497          premieryf = y1 
    498          dernierxf = x2 
    499          dernieryf = y2 
    500       ENDIF ELSE BEGIN 
    501          dom =   where( (glamf ge lon1) $ 
    502                         and (glamf le lon2) $ 
    503                         and (gphif ge lat1) $ 
    504                         and (gphif le lat2) ) 
    505          if (dom[0] eq -1) then begin 
    506             if keyword_set(findalways) then begin 
    507                domdef, prof1,prof2, GRILLE=grille, MEMEINDICES = memeindices, _extra = ex 
    508                return 
    509             ENDIF ELSE BEGIN 
    510                ras = report('ATTENTION, la boite ne contient pas de points f') 
    511                return 
    512             ENDELSE 
    513          ENDIF 
    514          ordonnee = dom/jpi 
    515          abscisse = dom-ordonnee*jpi 
    516          ordonnee = ordonnee(uniq(ordonnee, sort(ordonnee))) 
    517          abscisse = abscisse(uniq(abscisse, sort(abscisse))) 
    518          premierxf = abscisse[0] 
    519          premieryf = ordonnee[0] 
    520          dernierxf = abscisse[n_elements(abscisse) -1] 
    521          dernieryf = ordonnee[n_elements(ordonnee) -1] 
    522       ENDELSE 
    523       nxf = dernierxf-premierxf+1 
    524       nyf = dernieryf-premieryf+1 
    525       if (total(glamf[premierxf:dernierxf,premieryf:dernieryf] $ 
    526                 NE glamf[premierxf:dernierxf,0]#replicate(1, nyf)) eq 0 $ 
    527           AND total(gphif[premierxf:dernierxf,premieryf:dernieryf] $ 
    528                     NE replicate(1, nxf)#(gphif[0,premieryf:dernieryf])[*]) eq 0 ) $ 
    529        then key_irregular = 0 ELSE key_irregular = 1 
    530    endif 
    531 ;------------------------------------------------------------------- 
    532    if keyword_set(key_performance) THEN print, 'temps domdef', systime(1)-tempsun  
     752    ENDIF 
     753;------------------------------------------------------------------- 
     754; vertical domain defined with the Z index 
     755;------------------------------------------------------------------- 
     756  ENDIF ELSE BEGIN  
     757    CASE N_PARAMS() OF 
     758      2:fstz = min([x1, x2], max = lstz) 
     759      4:return 
     760      6:fstz = min([z1, z2], max = lstz) 
     761    ENDCASE 
     762    IF fstz LT 0 OR lstz GE jpk THEN BEGIN 
     763      ras = report('Bad definition of X1, X2, Z1 or Z2') 
     764      return 
     765    ENDIF 
     766    nz = lstz - fstz + 1 
     767; find vert1t, vert2t, firstzt, firstzt, nzt 
     768; according to (x1, x2) or (z1, z2) 
     769    IF (where(gridtype eq 'T'))[0] NE -1 THEN BEGIN 
     770      vert1t = min(gdept[fstz:lstz], max = vert2t) 
     771      firstzt = fstz & lastzt = lstz & nzt = nz 
     772    ENDIF       
     773; find vert1w, vert2w, firstzw, firstzw, nzw 
     774; according to (x1, x2) or (z1, z2) 
     775    IF (where(gridtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN BEGIN 
     776      vert1w = min(gdepw[fstz:lstz], max = vert2w) 
     777      firstzw = fstz & lastzw = lstz & nzw = nz 
     778    ENDIF       
     779    vert1 = min([vert1t, vert1w]) 
     780    vert2 = max([vert2t, vert2w]) 
     781  ENDELSE 
     782; 
     783;------------------------------------------------------------------- 
     784  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     785@updateold 
     786  ENDIF  
     787;------------------------------------------------------------------- 
     788  if keyword_set(key_performance) THEN print, 'temps domdef', systime(1)-tempsun  
    533789;------------------------------------------------------------------- 
    534790 
    535791;------------------------------------------------------------------- 
    536    return 
     792  return 
    537793end 

  • trunk/ToBeReviewed/GRILLE/f2v.pro

    r12 r13  
    3838;------------------------------------------------------------ 
    3939FUNCTION f2v, temp 
    40 @common 
     40;--------------------------------------------------------- 
     41@cm_4mesh 
     42@cm_4data 
     43@cm_4cal 
     44  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     45@updatenew 
     46  ENDIF 
     47;--------------------------------------------------------- 
    4148   res = temp 
    4249;on force nxt=nxf, etc ... 
    43    premierxv = premierxf 
    44    dernierxv = dernierxf 
    45    premieryv = premieryf 
    46    dernieryv = dernieryf 
     50   firstxv = firstxf 
     51   lastxv = lastxf 
     52   firstyv = firstyf 
     53   lastyv = lastyf 
    4754   nxv = nxf 
    4855   nyv = nyf 
    4956   vargrid = 'V' 
    5057   if NOT keyword_set(valmask) then valmask = 1e20 
    51    lon1 = glamv[premierxv, 0] 
    52    lon2 = glamf[dernierxf, 0] 
     58   lon1 = glamv[firstxv, 0] 
     59   lon2 = glamf[lastxf, 0] 
    5360 
    5461; cas sur la taille du tableau et application  
    … …  
    6067            taille[1] eq nxf and taille[2] eq nyf: 
    6168            taille[1] eq jpi and taille[2] eq jpj: $ 
    62              res=res[premierxf:dernierxf, premieryf:dernieryf] 
     69             res=res[firstxf:lastxf, firstyf:lastyf] 
    6370            else: $ 
    6471             return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 
    6572         endcase 
    66          mask = (fmask())[premierxf:dernierxf, premieryf:dernieryf, 0] 
     73         mask = (fmask())[firstxf:lastxf, firstyf:lastyf, 0] 
    6774         terre = where(mask EQ 0)  
    6875         IF terre[0] NE -1 THEN res[terre] = !values.f_nan 
    6976         res = 0.5*(res + shift(res, 1, 0)) 
    70          if NOT (keyword_set(key_periodique) AND nxf EQ jpi) then res[0, *] = !values.f_nan 
    71          mask = (vmask())[premierxf:dernierxf, premieryf:dernieryf, 0] 
     77         if NOT (keyword_set(key_periodic) AND nxf EQ jpi) then res[0, *] = !values.f_nan 
     78         mask = (vmask())[firstxf:lastxf, firstyf:lastyf, 0] 
    7279         terre = where(mask EQ 0)  
    7380         IF terre[0] NE -1 THEN res[terre] = valmask 
    … …  
    7784            taille[1] eq nxf and taille[2] eq nyf AND taille[3] EQ nzt: 
    7885            taille[1] eq nxf and taille[2] eq nyf AND taille[3] EQ jpk: $ 
    79              res=res[*, *, premierzt:dernierzt] 
     86             res=res[*, *, firstzt:lastzt] 
    8087            taille[1] eq nxf and taille[2] eq nyf AND taille[3] EQ jpt: 
    8188            taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk: $ 
    82              res=res[premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt] 
     89             res=res[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt] 
    8390            taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpt: $ 
    84              res=res[premierxf:dernierxf, premieryf:dernieryf, *] 
     91             res=res[firstxf:lastxf, firstyf:lastyf, *] 
    8592            else: $ 
    8693             return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 
    8794         ENDCASE 
    8895         if taille[3] EQ jpt then begin 
    89             mask = (fmask())[premierxf:dernierxf, premieryf:dernieryf, dernierzt*(nzt NE jpk)] 
     96            mask = (fmask())[firstxf:lastxf, firstyf:lastyf, lastzt*(nzt NE jpk)] 
    9097            mask = temporary(mask[*])#replicate(1, jpt) 
    9198            mask = reform(mask, nxf, nyf, jpt, /over) 
    92          ENDIF ELSE mask = (fmask())[premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt] 
     99         ENDIF ELSE mask = (fmask())[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt] 
    93100         terre = where(temporary(mask) EQ 0)  
    94101         IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan 
    95102         res = 0.5*(res + shift(res, 1, 0, 0)) 
    96          if NOT (keyword_set(key_periodique) AND nxf EQ jpi) then res[0, *, *] = !values.f_nan 
     103         if NOT (keyword_set(key_periodic) AND nxf EQ jpi) then res[0, *, *] = !values.f_nan 
    97104         if taille[3] EQ jpt then BEGIN 
    98             mask = tmask[premierxf:dernierxf, premieryf:dernieryf, dernierzt*(nzt NE jpk)] 
     105            mask = tmask[firstxf:lastxf, firstyf:lastyf, lastzt*(nzt NE jpk)] 
    99106            mask = temporary(mask[*])#replicate(1, jpt) 
    100107            mask = reform(mask, nxf, nyf, jpt, /over) 
    101          ENDIF ELSE mask = (vmask())[premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt] 
     108         ENDIF ELSE mask = (vmask())[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt] 
    102109         terre = where(temporary(mask) EQ 0)  
    103110         IF terre[0] NE -1 THEN res[temporary(terre)] = valmask 
    … …  
    107114            taille[1] eq nxf and taille[2] eq nyf AND taille[3] EQ nzt AND taille[4] EQ jpt: 
    108115            taille[1] eq nxf and taille[2] eq nyf AND taille[3] EQ jpk AND taille[4] EQ jpt: $ 
    109              res=res[*, *, premierzt:dernierzt, *] 
     116             res=res[*, *, firstzt:lastzt, *] 
    110117            taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk AND taille[4] EQ jpt: $ 
    111              res=res[premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt, *] 
     118             res=res[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt, *] 
    112119            else: $ 
    113120             return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 
    114121         ENDCASE 
    115          mask = (fmask())[premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt] 
     122         mask = (fmask())[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt] 
    116123         mask = temporary(mask[*])#replicate(1, jpt) 
    117124         mask = reform(mask, nxf, nyf, nzt, jpt, /over) 
    … …  
    119126         IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan 
    120127         res = 0.5*(res + shift(res, 1, 0, 0, 0)) 
    121          if NOT (keyword_set(key_periodique) AND nxf EQ jpi) then res[0, *, *, *] = !values.f_nan 
    122          mask = (vmask())[premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt] 
     128         if NOT (keyword_set(key_periodic) AND nxf EQ jpi) then res[0, *, *, *] = !values.f_nan 
     129         mask = (vmask())[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt] 
    123130         mask = temporary(mask[*])#replicate(1, jpt) 
    124131         mask = reform(mask, nxf, nyf, nzt, jpt, /over) 
    … …  
    128135   endcase 
    129136 
     137   IF NOT keyword_set(key_forgetold) THEN BEGIN 
     138   @updateold 
     139   ENDIF  
     140    
    130141   return, res 
    131142END 

  • trunk/ToBeReviewed/GRILLE/fmask.pro

    r12 r13  
    2525;------------------------------------------------------------ 
    2626FUNCTION fmask 
    27 @common 
    28    tempsun = systime(1)         ; pour key_performance 
    29    if jpk EQ 1 then begin 
    30       res = tmask*shift(tmask, -1, 0)*shift(tmask, 0, -1)*shift(tmask, -1, -1) 
    31       if NOT keyword_set(key_periodique) then res[jpi-1, *] = fmaskredy 
    32       res[*, jpj-1] = fmaskredx 
    33    ENDIF ELSE BEGIN 
    34       res = tmask*shift(tmask, -1, 0, 0)*shift(tmask, 0, -1, 0)*shift(tmask, -1, -1, 0) 
    35       if NOT keyword_set(key_periodique) then res[jpi-1, *, *] = fmaskredy 
    36       res[*, jpj-1, *] = fmaskredx 
    37    ENDELSE 
    38    if keyword_set(key_performance) THEN print, 'temps fmask', systime(1)-tempsun  
    39  
    40    return, res 
     27;--------------------------------------------------------- 
     28@cm_4mesh 
     29  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     30@updatenew 
     31  ENDIF 
     32;--------------------------------------------------------- 
     33  tempsun = systime(1)          ; pour key_performance 
     34;  
     35  CASE size(tmask, /n_dimensions) OF 
     36    2:res = tmask*shift(tmask, -1, 0)*shift(tmask, 0, -1)*shift(tmask, -1, -1) 
     37    3:res = tmask*shift(tmask, -1, 0, 0)*shift(tmask, 0, -1, 0)*shift(tmask, -1, -1, 0) 
     38  ENDCASE 
     39; 
     40  if NOT keyword_set(key_periodic) then res[jpi-1, *, *] = fmaskredy 
     41  res[*, jpj-1, *] = fmaskredx 
     42; 
     43  if keyword_set(key_performance) THEN print, 'temps fmask', systime(1)-tempsun  
     44   
     45  return, res 
    4146end 

  • trunk/ToBeReviewed/GRILLE/grille.pro

    r12 r13  
    1313; 
    1414; CALLING SEQUENCE: 
    15 ;  grille,mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz,e1,e2,e3 
     15;  grille,mask,glam,gphi,gdep,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz,e1,e2,e3 
    1616; 
    1717; INPUTS:rien. ATTENTION les choix de la grille se fait a partir de la 
    … …  
    3232;         /NOTRI: utile seulement qd TRI est active. dans ce cas 
    3333;         grille retourne -1 ds la variable tri meme si la variable du 
    34 ;         common triangles est definie et differente de -1 
    35 ; 
    36 ; OUTPUTS:mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz, 
    37 ;         dernierx,derniery,dernierz,e1,e2,e3 
     34;         common triangles_list est definie et differente de -1 
     35; 
     36;         /WDEPTH: to specify that the field is at W depth instad of T  
     37;         depth (automatically activated if vargrid eq 'W') 
     38; 
     39; OUTPUTS:mask,glam,gphi,gdep,nx,ny,nz,firstx,firsty,firstz, 
     40;         lastx,lasty,lastz,e1,e2,e3 
    3841; 
    3942;         pour leur definition cf domdef et la gestion des sous 
    … …  
    4346;         mask, glam et gphi il suffit de taper grille, mask, glam, gphi 
    4447; 
    45 ; COMMON BLOCKS: 
    46 ;       common.pro      congridseb.pro 
     48; COMMON BLOCKS: cm_4mesh and cm_4data 
    4749; 
    4850; SIDE EFFECTS: utilise la variable globale vargird 
    … …  
    5961;------------------------------------------------------------ 
    6062;------------------------------------------------------------ 
    61 pro grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz,dernierx, derniery, dernierz, e1,e2,e3, TRI = tri, NOTRI = notri, TOUT = tout, FORPLT = forplt, _EXTRA = ex 
    62 @common 
    63    tempsun = systime(1)         ; pour key_performance 
    64 ;------------------------------------------------------------ 
    65    if keyword_set(tout) then begin 
    66       oldboite = [lon1, lon2, lat1, lat2, prof1, prof2] 
    67       domdef, grille = vargrid, _EXTRA = ex 
    68    endif 
    69    tempdeux = systime(1)        ; pour key_performance =2 
    70    CASE 1 OF  
    71 ;------------------------------------------------------------ 
    72 ; grille T 
    73 ;------------------------------------------------------------ 
    74       vargrid eq 'OPAPTDHT' or vargrid eq 'OPAPT3DT' $ 
    75        or vargrid eq 'T': begin 
     63pro grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz, e1, e2, e3, TRI = tri, NOTRI = notri, TOUT = tout, FORPLT = forplt, IFPLTZ = ifpltz, WDEPTH = wdepth, _EXTRA = ex 
     64;------------------------------------------------------------ 
     65; include commons 
     66@cm_4mesh 
     67@cm_4data 
     68  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     69@updatenew 
     70  ENDIF 
     71;--------------------- 
     72  tempsun = systime(1)          ; pour key_performance 
     73;------------------------------------------------------------ 
     74  vargrid = strupcase(strmid(vargrid,0,/reverse_offset)) 
     75; 
     76  if vargrid eq 'W' then wdepth = 1 
     77  if keyword_set(tout) then begin 
     78    savedbox = 1b 
     79    saveboxparam, 'boxparam4grille.dat' 
     80    domdef, gridtype = vargrid, _EXTRA = ex 
     81  endif 
     82  tempdeux = systime(1)         ; pour key_performance =2 
     83;------------------------------------------------------------ 
     84;------------------------------------------------------------ 
     85  IF keyword_set(wdepth) THEN BEGIN 
     86    firstz = firstzw  
     87    lastz = lastzw 
     88    nz = nzw 
     89  ENDIF ELSE BEGIN 
     90    firstz = firstzt 
     91    lastz = lastzt 
     92    nz = nzt 
     93  ENDELSE 
     94;------------------------------------------------------------ 
     95;------------------------------------------------------------ 
     96  CASE 1 OF  
     97;------------------------------------------------------------ 
     98; grille T and W 
     99;------------------------------------------------------------ 
     100    vargrid eq 'T' OR vargrid eq 'W' : begin 
    76101;scalaires 
    77          nx=nxt 
    78          ny=nyt 
    79          nz=nzt 
    80          premierx = premierxt 
    81          premiery = premieryt 
    82          premierz = premierzt 
    83          dernierx = dernierxt 
    84          derniery = dernieryt 
    85          dernierz = dernierzt 
     102      nx = nxt 
     103      ny = nyt 
     104      firstx = firstxt 
     105      firsty = firstyt 
     106      lastx = lastxt 
     107      lasty = lastyt 
    86108;vecteurs 2d 
    87          glam=glamt[premierx:dernierx, premiery:derniery] 
    88          gphi=gphit[premierx:dernierx, premiery:derniery] 
    89          e1  =e1t[premierx:dernierx, premiery:derniery] 
    90          e2  =e2t[premierx:dernierx, premiery:derniery] 
     109      IF n_elements(glam) NE 1 THEN glam = glamt[firstx:lastx, firsty:lasty] 
     110      IF n_elements(gphi) NE 1 THEN gphi = gphit[firstx:lastx, firsty:lasty] 
     111      IF n_elements(e1) NE 1 THEN e1 = e1t[firstx:lastx, firsty:lasty] 
     112      IF n_elements(e2) NE 1 THEN e2 = e2t[firstx:lastx, firsty:lasty] 
    91113;vecteurs 3d 
    92          mask=tmask[premierx:dernierx, premiery:derniery, premierz:dernierz] 
    93       end 
    94 ;------------------------------------------------------------ 
    95 ; grille W 
    96 ;------------------------------------------------------------ 
    97       vargrid eq 'OPAPT3DW' or vargrid eq 'W': begin 
     114      IF keyword_set(forplt) THEN mask = tmask[firstx:lastx, firsty:lasty, firstz] $ 
     115      ELSE IF n_elements(mask) NE 1 THEN mask = tmask[firstx:lastx, firsty:lasty, firstz:lastz] 
     116    end 
     117;------------------------------------------------------------ 
     118; grille U 
     119;------------------------------------------------------------ 
     120    vargrid eq 'U': begin 
    98121;scalaires 
    99          nx=nxt 
    100          ny=nyt 
    101          nz=nzw 
    102          premierx = premierxt 
    103          premiery = premieryt 
    104          premierz = premierzw 
    105          dernierx = dernierxt 
    106          derniery = dernieryt 
    107          dernierz = dernierzw 
     122      nx = nxu 
     123      ny = nyu 
     124      firstx = firstxu 
     125      firsty = firstyu 
     126      lastx = lastxu 
     127      lasty = lastyu 
    108128;vecteurs 2d 
    109          terre = where(tmask[*, *, 0] EQ 0) 
    110          glam=glamt[premierx:dernierx, premiery:derniery] 
    111          gphi=gphit[premierx:dernierx, premiery:derniery] 
    112          e1  =e1t[premierx:dernierx, premiery:derniery] 
    113          e2  =e2t[premierx:dernierx, premiery:derniery] 
     129      IF n_elements(glam) NE 1 THEN glam = glamu[firstx:lastx, firsty:lasty] 
     130      IF n_elements(gphi) NE 1 THEN gphi = gphiu[firstx:lastx, firsty:lasty] 
     131      if keyword_set(forplt) then BEGIN 
     132        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz] 
     133        eastboarder = mask-shift(mask, 1, 0)*mask 
     134        westboarder = mask-shift(mask, -1, 0)*mask 
     135        if key_periodic NE 1 OR nx NE jpi then westboarder[nx-1, *] = 0b 
     136        tmp1 = shift(eastboarder, 0, 1) 
     137        tmp1[*, 0] = 0b 
     138        tmp2 = shift(eastboarder, 0, -1) 
     139        tmp2[*, ny-1] = 0b 
     140        add = (temporary(tmp1)+temporary(tmp2))*(1b-eastboarder)*(1b-temporary(westboarder)) 
     141        eastboarder = temporary(eastboarder)+temporary(add) 
     142        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b 
     143        tmp1[*, ny-1] = 1b 
     144        tmp1[*, 0] = 1b 
     145        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b 
     146        if key_periodic NE 1 OR nx NE jpi then begin 
     147          tmp2[nx-1, *] = 1b 
     148          tmp2[0, *] = 0b 
     149        endif 
     150        no1 = temporary(tmp1)*temporary(tmp2) 
     151        tmp = temporary(eastboarder)*temporary(no1)*mask 
     152        mask[0:nx-2, *] = 0b  
     153        tmp = temporary(tmp)+temporary(mask) 
     154        tmp = where(tmp GE 1) 
     155        if tmp[0] NE -1 then begin 
     156          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp] 
     157          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp] 
     158        endif 
     159      ENDIF 
     160      IF n_elements(e1) NE 1 THEN e1  = e1u[firstx:lastx, firsty:lasty] 
     161      IF n_elements(e2) NE 1 THEN e2  = e2u[firstx:lastx, firsty:lasty] 
    114162;vecteurs 3d 
    115          mask=tmask[premierx:dernierx, premiery:derniery, premierz:dernierz] 
    116       end 
    117 ;------------------------------------------------------------ 
    118 ; grille U 
    119 ;------------------------------------------------------------ 
    120       vargrid eq 'OPAPTDHU' or vargrid eq 'OPAPT3DU' $ 
    121        or vargrid eq 'U': begin 
     163      IF keyword_set(forplt) THEN mask = (umask())[firstx:lastx, firsty:lasty, firstz] $ 
     164      ELSE IF n_elements(mask) NE 1 THEN mask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz] 
     165    end 
     166;------------------------------------------------------------ 
     167; grille V 
     168;------------------------------------------------------------ 
     169    vargrid eq 'OPAPTDHV' or vargrid eq 'OPAPT3DV' $ 
     170      or vargrid eq 'V': begin 
    122171;scalaires 
    123          nx=nxu 
    124          ny=nyu 
    125          nz=nzt 
    126          premierx = premierxu 
    127          premiery = premieryu 
    128          premierz = premierzt 
    129          dernierx = dernierxu 
    130          derniery = dernieryu 
    131          dernierz = dernierzt 
     172      nx = nxv 
     173      ny = nyv 
     174      firstx = firstxv 
     175      firsty = firstyv 
     176      lastx = lastxv 
     177      lasty = lastyv 
    132178;vecteurs 2d 
    133          glam=glamu[premierx:dernierx, premiery:derniery] 
    134          gphi=gphiu[premierx:dernierx, premiery:derniery] 
    135          if keyword_set(forplt) then BEGIN  
    136             mask = tmask[premierx:dernierx, premiery:derniery, niveau-1] 
    137             terre = where(mask EQ 0) 
    138             if terre[0] NE -1 then begin 
    139                glam[terre] = (glamt[premierx:dernierx, premiery:derniery])[terre] 
    140                gphi[terre] = (gphit[premierx:dernierx, premiery:derniery])[terre] 
    141             endif 
    142          ENDIF 
    143          e1  =e1u[premierx:dernierx, premiery:derniery] 
    144          e2  =e2u[premierx:dernierx, premiery:derniery] 
     179      IF n_elements(glam) NE 1 THEN glam = glamv[firstx:lastx, firsty:lasty] 
     180      IF n_elements(gphi) NE 1 THEN gphi = gphiv[firstx:lastx, firsty:lasty] 
     181      if keyword_set(forplt) then BEGIN  
     182        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz] 
     183        northboarder = mask-shift(mask, 0, 1)*mask 
     184        southboarder = mask-shift(mask, 0, -1)*mask 
     185        southboarder[*, ny-1] = 0b 
     186        tmp1 = shift(northboarder, -1, 0) 
     187        if key_periodic NE 1 OR nx NE jpi then tmp1[nx-1, *] = 0b 
     188        tmp2 = shift(northboarder, 1, 0) 
     189        if key_periodic NE 1 OR nx NE jpi then tmp2[0, *] = 0b 
     190        add = (temporary(tmp1)+temporary(tmp2))*(1b-northboarder)*(1b-southboarder) 
     191        northboarder = temporary(northboarder)+temporary(add) 
     192        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b 
     193        tmp1[*, ny-1] = 1b 
     194        tmp1[*, 0] = 0b 
     195        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b 
     196        if key_periodic NE 1 OR nx NE jpi then begin 
     197          tmp2[nx-1, *] = 1b 
     198          tmp2[0, *] = 1b 
     199        endif 
     200        no1 = temporary(tmp1)*temporary(tmp2) 
     201        tmp = temporary(northboarder)*mask*temporary(no1) 
     202        mask[*, 0:ny-2] = 0b 
     203        tmp = temporary(tmp)+temporary(mask) 
     204        tmp = where(tmp GE 1) 
     205        if tmp[0] NE -1 then begin 
     206          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp] 
     207          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp] 
     208        endif 
     209      ENDIF 
     210      IF n_elements(e1) NE 1 THEN e1  = e1v[firstx:lastx, firsty:lasty] 
     211      IF n_elements(e2) NE 1 THEN e2  = e2v[firstx:lastx, firsty:lasty] 
    145212;vecteurs 3d 
    146          mask=(umask())[premierx:dernierx, premiery:derniery, premierz:dernierz] 
    147       end 
    148 ;------------------------------------------------------------ 
    149 ; grille V 
    150 ;------------------------------------------------------------ 
    151       vargrid eq 'OPAPTDHV' or vargrid eq 'OPAPT3DV' $ 
    152        or vargrid eq 'V': begin 
     213      IF keyword_set(forplt) THEN mask = (vmask())[firstx:lastx, firsty:lasty, firstz] $ 
     214      ELSE IF n_elements(mask) NE 1 THEN mask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz] 
     215    end 
     216;------------------------------------------------------------ 
     217; grille F 
     218;------------------------------------------------------------ 
     219    vargrid eq 'OPAPTDHF' or vargrid eq 'OPAPT3DF' $  
     220      or vargrid eq 'F': begin 
    153221;scalaires 
    154          nx=nxv 
    155          ny=nyv 
    156          nz=nzt 
    157          premierx = premierxv 
    158          premiery = premieryv 
    159          premierz = premierzt 
    160          dernierx = dernierxv 
    161          derniery = dernieryv 
    162          dernierz = dernierzt 
     222      nx = nxf 
     223      ny = nyf 
     224      firstx = firstxf 
     225      firsty = firstyf 
     226      lastx = lastxf 
     227      lasty = lastyf 
    163228;vecteurs 2d 
    164          glam=glamv[premierx:dernierx, premiery:derniery] 
    165          gphi=gphiv[premierx:dernierx, premiery:derniery] 
    166          if keyword_set(forplt) then BEGIN  
    167             mask = tmask[premierx:dernierx, premiery:derniery, niveau-1] 
    168             terre = where(mask EQ 0) 
    169             if terre[0] NE -1 then begin 
    170                glam[terre] = (glamt[premierx:dernierx, premiery:derniery])[terre] 
    171                gphi[terre] = (gphit[premierx:dernierx, premiery:derniery])[terre] 
    172             endif 
    173          ENDIF 
    174          e1  =e1v[premierx:dernierx, premiery:derniery] 
    175          e2  =e2v[premierx:dernierx, premiery:derniery] 
     229      IF n_elements(glam) NE 1 THEN glam = glamf[firstx:lastx, firsty:lasty] 
     230      IF n_elements(gphi) NE 1 THEN gphi = gphif[firstx:lastx, firsty:lasty] 
     231      if keyword_set(forplt) then BEGIN  
     232        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz] 
     233        eastboarder = mask-shift(mask, 1, 0)*mask 
     234        westboarder = mask-shift(mask, -1, 0)*mask 
     235        westboarder[nx-1, *] = 0b 
     236        northboarder = mask-shift(mask, 0, 1)*mask 
     237        southboarder = mask-shift(mask, 0, -1)*mask 
     238        southboarder[*, ny-1] = 0b 
     239        tmp1 = shift(northboarder, -1, 0) 
     240        if key_periodic NE 1 OR nx NE jpi then tmp1[nx-1, *] = 0b 
     241        tmp2 = shift(northboarder, 1, 0) 
     242        if key_periodic NE 1 OR nx NE jpi then tmp2[0, *] = 0b 
     243        add = (temporary(tmp1)+temporary(tmp2))*(1b-northboarder)*(1b-southboarder) 
     244        northboarder = temporary(northboarder)+temporary(add) 
     245        tmp1 = shift(eastboarder, 0, 1) 
     246        tmp1[*, 0] = 0b 
     247        tmp2 = shift(eastboarder, 0, -1) 
     248        tmp2[*, ny-1] = 0b 
     249        add = (temporary(tmp1)+temporary(tmp2))*(1b-eastboarder)*(1b-temporary(westboarder)) 
     250        eastboarder = temporary(eastboarder)+temporary(add) 
     251        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b 
     252        tmp1[*, ny-1] = 1b 
     253        tmp1[*, 0] = 1b 
     254        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b 
     255        if key_periodic NE 1 OR nx NE jpi then begin 
     256          tmp2[nx-1, *] = 1b 
     257          tmp2[0, *] = 1b 
     258        endif 
     259        no1 = temporary(tmp1)*temporary(tmp2) 
     260        tmp = (temporary(northboarder)+temporary(eastboarder))*mask*temporary(no1) 
     261        mask[0:nx-2, *] = 0b  
     262        mask[*, 0:ny-2] = 0b  
     263        tmp = temporary(tmp)+temporary(mask) 
     264        tmp = where(tmp GE 1) 
     265        if tmp[0] NE -1 then begin 
     266          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp] 
     267          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp] 
     268        endif 
     269      ENDIF 
     270      IF n_elements(e1) NE 1 THEN e1  = e1f[firstx:lastx, firsty:lasty] 
     271      IF n_elements(e2) NE 1 THEN e2  = e2f[firstx:lastx, firsty:lasty] 
    176272;vecteurs 3d 
    177          mask=(vmask())[premierx:dernierx, premiery:derniery, premierz:dernierz] 
    178       end 
    179 ;------------------------------------------------------------ 
    180 ; grille F 
    181 ;------------------------------------------------------------ 
    182       vargrid eq 'OPAPTDHF' or vargrid eq 'OPAPT3DF' $  
    183        or vargrid eq 'F': begin 
    184 ;scalaires 
    185          nx=nxf 
    186          ny=nyf 
    187          nz=nzt 
    188          premierx = premierxf 
    189          premiery = premieryf 
    190          premierz = premierzt 
    191          dernierx = dernierxf 
    192          derniery = dernieryf 
    193          dernierz = dernierzt 
    194 ;vecteurs 2d 
    195          glam=glamf[premierx:dernierx, premiery:derniery] 
    196          gphi=gphif[premierx:dernierx, premiery:derniery] 
    197          if keyword_set(forplt) then BEGIN  
    198             mask = tmask[premierx:dernierx, premiery:derniery, niveau-1] 
    199             terre = where(mask EQ 0) 
    200             if terre[0] NE -1 then begin 
    201                glam[terre] = (glamt[premierx:dernierx, premiery:derniery])[terre] 
    202                gphi[terre] = (gphit[premierx:dernierx, premiery:derniery])[terre] 
    203             endif 
    204          ENDIF 
    205          e1  =e1f[premierx:dernierx, premiery:derniery] 
    206          e2  =e2f[premierx:dernierx, premiery:derniery] 
    207 ;vecteurs 3d 
    208          mask=(fmask())[premierx:dernierx, premiery:derniery, premierz:dernierz] 
     273      IF keyword_set(forplt) THEN mask = (fmask())[firstx:lastx, firsty:lasty, firstz] $ 
     274      ELSE IF n_elements(mask) NE 1 THEN mask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz] 
     275    END 
     276;------------------------------------------------------------ 
     277    ELSE:BEGIN  
     278      ras = report('Wrong definition of vargrid = '+vargrid+'. Only T, U, V, W or F are acceptable') 
     279      stop 
     280    END 
     281  ENDCASE 
     282  IF testvar(var = key_performance) EQ 2 THEN $ 
     283    print, 'temps grille: attribution des scalaires, vecteurs et tableaux ', systime(1)-tempdeux 
     284; 
     285;------------------------------------------------------------ 
     286;------------------------------------------------------------ 
     287;------------------------------------------------------------ 
     288; Variables se rapportant a la dimension verticale 
     289;------------------------------------------------------------ 
     290;------------------------------------------------------------ 
     291;------------------------------------------------------------ 
     292; 
     293; 
     294  tempdeux = systime(1)         ; pour key_performance =2 
     295  if keyword_set(wdepth) then begin 
     296    gdep = gdepw[firstz:lastz] 
     297    e3 = e3w[firstz:lastz] 
     298  endif else begin 
     299    gdep = gdept[firstz:lastz] 
     300    e3 = e3t[firstz:lastz] 
     301  ENDELSE 
     302; for the vertical sections with partial steps 
     303  IF keyword_set(ifpltz) AND keyword_set(key_partialstep) THEN BEGIN 
     304    CASE 1 OF 
     305      ifpltz EQ 'xz' AND ny EQ 1:BEGIN 
     306        bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3) 
     307        good = where(bottom NE 0 AND bottom NE nz+keyword_set(wdepth)) 
     308        bottom = lindgen(nx)+(bottom-1l+keyword_set(wdepth))*nx 
     309        IF good[0] NE -1 THEN BEGIN 
     310          bottom = bottom[good] 
     311          IF lastz EQ jpk-1 THEN gdep[nz-1] = max(hdepw) 
     312          gdep = replicate(1, nx)#gdep 
     313          if keyword_set(wdepth) THEN $ 
     314            truegdep = hdepw[firstx:lastx, firsty:lasty] $ 
     315          ELSE truegdep = hdept[firstx:lastx, firsty:lasty] 
     316          gdep[bottom] = truegdep[good] 
     317        ENDIF 
    209318      END 
    210 ;------------------------------------------------------------ 
    211       ELSE:BEGIN  
    212          ras = report('Vauvaise definition de la variable vargrid: '+vargrid+'ceete variable doit etre T, U, V, W ou F') 
    213          stop 
     319      ifpltz EQ 'yz' AND nx EQ 1:BEGIN 
     320        bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3) 
     321        good = where(bottom NE 0 AND bottom NE nz+keyword_set(wdepth)) 
     322        bottom = lindgen(ny)+(bottom-1l+keyword_set(wdepth))*ny 
     323        IF good[0] NE -1 THEN BEGIN 
     324          bottom = bottom[good] 
     325          IF lastz EQ jpk-1 THEN gdep[nz-1] = max(hdepw) 
     326          gdep = replicate(1, ny)#gdep 
     327          if keyword_set(wdepth) THEN $ 
     328            truegdep = hdepw[firstx:lastx, firsty:lasty] $ 
     329          ELSE truegdep = hdept[firstx:lastx, firsty:lasty] 
     330          gdep[bottom] = truegdep[good] 
     331        ENDIF 
    214332      END 
    215    ENDCASE 
    216    IF testvar(var = key_performance) EQ 2 THEN $ 
    217     print, 'temps grille: attribution des scalaires, vecteurs et tableaux ', systime(1)-tempdeux 
    218 ;------------------------------------------------------------ 
    219 ; Variables se rapportant a la dimension verticale 
    220 ;------------------------------------------------------------ 
    221    tempdeux = systime(1)        ; pour key_performance =2 
    222    if vargrid eq 'OPAPT3DW' or vargrid eq 'W' then begin 
    223       gdep = gdepw[premierz:dernierz] 
    224       e3=e3w[premierz:dernierz] 
    225    endif else begin 
    226       gdep = gdept[premierz:dernierz] 
    227       e3=e3t[premierz:dernierz] 
    228    endelse 
    229    IF testvar(var = key_performance) EQ 2 THEN $ 
     333      ELSE: 
     334    ENDCASE 
     335  ENDIF 
     336  IF testvar(var = key_performance) EQ 2 THEN $ 
    230337    print, 'temps grille: Variables se rapportant a la dimension verticale ', systime(1)-tempdeux 
    231338;------------------------------------------------------------ 
    232339; vecteur triangulation Qd TRI est active 
    233340;------------------------------------------------------------ 
    234    if arg_present(TRI) then $ 
    235     if triangles[0] EQ -1 OR keyword_set(notri) then tri = -1 ELSE BEGIN  
    236       tempdeux = systime(1)     ; pour key_performance =2 
    237       msk = bytarr(jpi, jpj) 
    238       msk[premierx:dernierx,premiery:derniery] = 1 
    239       ind = where( msk[triangles[0, *]]*msk[triangles[1, *]]*msk[triangles[2, *]] EQ 1 ) 
    240       tri =triangles[*, ind]-(premierx+premiery*jpi) 
    241       y = tri/jpi 
    242       x = tri-y*jpi 
    243       tri = x+y*nx 
    244       IF testvar(var = key_performance) EQ 2 THEN $ 
    245        print, 'temps grille: decoupage de la triangulation ', systime(1)-tempdeux 
    246    ENDELSE 
     341  if arg_present(TRI) then $ 
     342    if triangles_list[0] EQ -1 OR keyword_set(notri) then tri = -1 ELSE BEGIN  
     343    tempdeux = systime(1)       ; pour key_performance =2 
     344    msk = bytarr(jpi, jpj) 
     345    msk[firstx:lastx, firsty:lasty] = 1 
     346    ind = where( msk[triangles_list[0, *]]*msk[triangles_list[1, *]]*msk[triangles_list[2, *]] EQ 1 ) 
     347    tri = triangles_list[*, ind]-(firstx+firsty*jpi) 
     348    y = tri/jpi 
     349    x = tri-y*jpi 
     350    tri = x+y*nx 
     351    IF testvar(var = key_performance) EQ 2 THEN $ 
     352      print, 'temps grille: decoupage de la triangulation ', systime(1)-tempdeux 
     353  ENDELSE 
    247354;------------------------------------------------------------------ 
    248355; pour s'assurer qu'il n'y a pas de dimension degenerees (=1) 
    … …  
    256363;    e3=reform(e3, /over) 
    257364 
    258    if keyword_set(tout) then domdef, oldboite, grille = vargrid 
    259    if keyword_set(key_performance) THEN print, 'temps grille', systime(1)-tempsun  
    260  
    261    return 
     365  if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grille.dat' 
     366  if keyword_set(key_performance) THEN print, 'temps grille', systime(1)-tempsun  
     367 
     368;------------------------------------------------------------ 
     369  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     370@updateold 
     371  ENDIF  
     372;--------------------- 
     373  return 
    262374 
    263375end 

  • trunk/ToBeReviewed/GRILLE/t2v.pro

    r12 r13  
    3838;------------------------------------------------------------ 
    3939FUNCTION t2v, temp 
    40 @common 
     40;--------------------------------------------------------- 
     41@cm_4mesh 
     42@cm_4data 
     43@cm_4cal 
     44  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     45@updatenew 
     46  ENDIF 
     47;--------------------------------------------------------- 
    4148   res = temp 
    4249 
    4350;on force nxt=nxv, etc ... 
    44    premierxv = premierxt 
    45    dernierxv = dernierxt 
    46    premieryv = premieryt 
    47    dernieryv = dernieryt 
     51   firstxv = firstxt 
     52   lastxv = lastxt 
     53   firstyv = firstyt 
     54   lastyv = lastyt 
    4855   nxv = nxt 
    4956   nyv = nyt 
    5057   vargrid = 'V' 
    5158   if NOT keyword_set(valmask) then valmask = 1e20 
    52    lat1 = gphit[0, premieryt] 
    53    lat2 = gphiv[0, dernieryv] 
     59   lat1 = gphit[0, firstyt] 
     60   lat2 = gphiv[0, lastyv] 
    5461 
    5562; cas sur la taille du tableau et application  
    … …  
    6168            taille[1] eq nxt and taille[2] eq nyt: 
    6269            taille[1] eq jpi and taille[2] eq jpj: $ 
    63              res=res[premierxt:dernierxt, premieryt:dernieryt] 
     70             res=res[firstxt:lastxt, firstyt:lastyt] 
    6471            else: $ 
    6572             return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 
    6673         endcase 
    67          mask = tmask[premierxt:dernierxt, premieryt:dernieryt, 0] 
     74         mask = tmask[firstxt:lastxt, firstyt:lastyt, 0] 
    6875         terre = where(mask EQ 0)  
    6976         IF terre[0] NE -1 THEN res[terre] = !values.f_nan 
    7077         res = 0.5*(res + shift(res, 0, -1)) 
    7178         res[*, nyt-1] = !values.f_nan 
    72          mask = (vmask())[premierxt:dernierxt, premieryt:dernieryt, 0] 
     79         mask = (vmask())[firstxt:lastxt, firstyt:lastyt, 0] 
    7380         terre = where(mask EQ 0)  
    7481         IF terre[0] NE -1 THEN res[terre] = valmask 
    … …  
    7885            taille[1] eq nxt and taille[2] eq nyt AND taille[3] EQ nzt: 
    7986            taille[1] eq nxt and taille[2] eq nyt AND taille[3] EQ jpk: $ 
    80              res=res[*, *, premierzt:dernierzt] 
     87             res=res[*, *, firstzt:lastzt] 
    8188            taille[1] eq nxt and taille[2] eq nyt AND taille[3] EQ jpt: 
    8289            taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk: $ 
    83              res=res[premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt] 
     90             res=res[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt] 
    8491            taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpt: $ 
    85              res=res[premierxt:dernierxt, premieryt:dernieryt, *] 
     92             res=res[firstxt:lastxt, firstyt:lastyt, *] 
    8693            else: $ 
    8794             return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 
    8895         ENDCASE 
    8996         if taille[3] EQ jpt then begin 
    90             mask = tmask[premierxt:dernierxt, premieryt:dernieryt, dernierzt*(nzt NE jpk)] 
     97            mask = tmask[firstxt:lastxt, firstyt:lastyt, lastzt*(nzt NE jpk)] 
    9198            mask = temporary(mask[*])#replicate(1, jpt) 
    9299            mask = reform(mask, nxt, nyt, jpt, /over) 
    93          ENDIF ELSE mask = tmask[premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt] 
     100         ENDIF ELSE mask = tmask[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt] 
    94101         terre = where(temporary(mask) EQ 0)  
    95102         IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan 
    … …  
    97104         res[*, nyt-1, *] = !values.f_nan 
    98105         if taille[3] EQ jpt then BEGIN 
    99             mask = (vmask())[premierxt:dernierxt, premieryt:dernieryt, dernierzt*(nzt NE jpk)] 
     106            mask = (vmask())[firstxt:lastxt, firstyt:lastyt, lastzt*(nzt NE jpk)] 
    100107            mask = temporary(mask[*])#replicate(1, jpt) 
    101108            mask = reform(mask, nxt, nyt, jpt, /over) 
    102          ENDIF ELSE mask = (vmask())[premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt] 
     109         ENDIF ELSE mask = (vmask())[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt] 
    103110         terre = where(temporary(mask) EQ 0)  
    104111         IF terre[0] NE -1 THEN res[temporary(terre)] = valmask 
    … …  
    108115            taille[1] eq nxt and taille[2] eq nyt AND taille[3] EQ nzt AND taille[4] EQ jpt: 
    109116            taille[1] eq nxt and taille[2] eq nyt AND taille[3] EQ jpk AND taille[4] EQ jpt: $ 
    110              res=res[*, *, premierzt:dernierzt, *] 
     117             res=res[*, *, firstzt:lastzt, *] 
    111118            taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk AND taille[4] EQ jpt: $ 
    112              res=res[premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt, *] 
     119             res=res[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt, *] 
    113120            else: $ 
    114121             return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 
    115122         ENDCASE 
    116          mask = tmask[premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt] 
     123         mask = tmask[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt] 
    117124         mask = temporary(mask[*])#replicate(1, jpt) 
    118125         mask = reform(mask, nxt, nyt, nzt, jpt, /over) 
    … …  
    121128         res = 0.5*(res + shift(res, 0, -1, 0, 0)) 
    122129         res[*, nyt-1, *, *] = !values.f_nan 
    123          mask = (vmask())[premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt] 
     130         mask = (vmask())[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt] 
    124131         mask = temporary(mask[*])#replicate(1, jpt) 
    125132         mask = reform(mask, nxt, nyt, nzt, jpt, /over) 
    … …  
    127134         IF terre[0] NE -1 THEN res[temporary(terre)] = valmask 
    128135      END 
    129    endcase 
     136    ENDCASE 
     137 
     138  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     139   @updateold 
     140  ENDIF  
    130141 
    131142   return, res 

  • trunk/ToBeReviewed/GRILLE/tracegrille.pro

    r12 r13  
    1111; CALLING SEQUENCE:tracegrille 
    1212;  
    13 ; INPUTS:none 
     13; INPUTS:glam et gphi, les tableaux 1d ou 2d des position en 
     14; longitude/latitude des points de la grille. Si glam et gphi ne sont 
     15; pas specifies, trace la grille specifiee par vargrid, sur le domaine 
     16; definit par le dernier domdef. 
    1417; 
    1518; KEYWORD PARAMETERS: 
    … …  
    2124;        qu''une ligne de j constant tout les ystride points 
    2225; 
    23 ;        OCEAN: pour ne tracer la grille que sur les points oceans 
     26;        /OCEAN: pour ne tracer la grille que sur les points oceans 
    2427;  
    25 ;        EARTH: pour ne tracer la grille que sur les points terre 
     28;        /EARTH: pour ne tracer la grille que sur les points terre 
     29; 
     30;        /RMOUT:select to remove all cell having one corner out of the 
     31;        plot boundaries (!x.range, !y.range) 
    2632; 
    2733;        + tous les mots clefs de la procedure PLOTS 
    … …  
    3137; COMMON BLOCKS:common.pro 
    3238; 
    33 ; SIDE EFFECTS:trace la grille specifiee par vargrid, sur le domaine 
    34 ; definit par le dernier domdef. 
     39; SIDE EFFECTS: 
    3540; 
    3641; RESTRICTIONS: 
    … …  
    3843; EXAMPLE: 
    3944; 
    40 ;     IDL> plt,indgen(jpi,jpj),/nocontour,/nocouleur 
     45;     IDL> plt,indgen(jpi,jpj),/nocontour,/nofill 
     46;     IDL> vargrid='T' 
    4147;     IDL> tracegrille,/ocean,color=20 
    4248;     IDL> tracegrille,/earth,color=80 
    … …  
    4955;------------------------------------------------------------ 
    5056;------------------------------------------------------------ 
    51 PRO tracegrille, OCEAN = ocean, EARTH = earth, XSTRIDE = xstride, YSTRIDE = ystride, _extra = extra 
    52 @common 
    53    tempsun = systime(1)         ; pour key_performance 
    54    if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c' 
     57PRO tracegrille, glamin, gphiin, OCEAN = ocean, EARTH = earth $ 
     58                 , XSTRIDE = xstride, YSTRIDE = ystride, RMOUT = rmout $ 
     59                 , _extra = extra 
     60;--------------------------------------------------------- 
     61@cm_4mesh 
     62@cm_4data 
     63  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     64@updatenew 
     65  ENDIF 
     66;--------------------------------------------------------- 
     67  tempsun = systime(1)          ; pour key_performance 
     68; to avoid warning message 
     69  oldexcept = !except 
     70  !except = 0 
     71  if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c' 
     72; 
     73  if n_elements(glamin) * n_elements(gphiin) EQ 0 then BEGIN 
     74    grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz 
     75    IF keyword_set(ocean) AND key_gridtype EQ 'c' THEN BEGIN 
     76; we reduce the mask to take into account the point located ON the coastline. 
     77      CASE vargrid OF 
     78        'U':BEGIN 
     79          mask = tmask[firstx:lastx, firsty:lasty] 
     80          IF NOT keyword_set(key_periodic) OR nx NE jpi $ 
     81            THEN tmpx = mask[nx-1, *] 
     82          mask = (mask+shift(mask, -1, 0)) < 1 
     83          IF NOT keyword_set(key_periodic) OR nx NE jpi $ 
     84            THEN mask[nx-1, *] = temporary(tmpx) 
     85        END 
     86        'V':BEGIN 
     87          mask = tmask[firstx:lastx, firsty:lasty] 
     88          tmpy = mask[*, ny-1] 
     89          mask = (mask+shift(mask, 0, -1)) < 1 
     90          mask[*, ny-1] = temporary(tmpy) 
     91        END 
     92        'F':BEGIN 
     93          mask = tmask[firstx:lastx, firsty:lasty] 
     94          IF NOT keyword_set(key_periodic) OR nx NE jpi $ 
     95            THEN tmpx = mask[nx-1, *] 
     96          tmpy = mask[*, ny-1] 
     97          mask = (mask+shift(mask, -1, 0)+shift(mask, 0, -1)+shift(mask, -1, -1)) < 1 
     98          mask[*, ny-1] = temporary(tmpy) 
     99          IF NOT keyword_set(key_periodic) OR nx NE jpi $ 
     100            THEN mask[nx-1, *] = temporary(tmpx) 
     101        END 
     102        ELSE: 
     103      ENDCASE 
     104    ENDIF 
     105  ENDIF ELSE BEGIN 
     106    glam = glamin 
     107    gphi = gphiin 
     108    IF (size(glam))[0] EQ 1 AND (size(gphi))[0] EQ 1 THEN BEGIN 
     109      nx = n_elements(glam)  
     110      ny = n_elements(gphi)  
     111      glam = glam#replicate(1, ny) 
     112      gphi = replicate(1, nx)#gphi 
     113    ENDIF ELSE BEGIN 
     114      nx = (size(glam))[1] 
     115      ny = (size(glam))[2] 
     116    ENDELSE 
     117  ENDELSE 
     118  if n_elements(mask) EQ 0 then mask = replicate(1b, nx, ny) 
     119  if (size(mask))[0] EQ 3 then mask = mask[*, *, 0] 
     120; 
     121  IF keyword_set(RMOUT) THEN BEGIN 
     122    out = where(glam GT max(!x.range) OR glam LT min(!x.range) $ 
     123                OR gphi GT max(!y.range) OR gphi LT min(!y.range)) 
     124    IF out[0] NE -1 THEN BEGIN 
     125      glam[out] = !values.f_nan 
     126      gphi[out] = !values.f_nan 
     127    ENDIF 
     128  ENDIF 
     129; 
     130  IF keyword_set(ocean) then BEGIN 
     131    earth = where(mask EQ 0) 
     132    if earth[0] NE -1 then begin 
     133      glam[earth] = !values.f_nan 
     134      gphi[earth] = !values.f_nan 
     135    ENDIF 
     136    earth = 0 
     137  ENDIF 
     138; 
     139  IF keyword_set(earth) THEN BEGIN 
     140    ocean = where(mask EQ 1) 
     141    if ocean[0] NE -1 then begin 
     142      glam[ocean] = !values.f_nan 
     143      gphi[ocean] = !values.f_nan 
     144    ENDIF 
     145    ocean = 0 
     146  ENDIF 
     147; 
     148  if NOT keyword_set(xstride) then xstride = 1 
     149  if NOT keyword_set(ystride) then ystride = 1 
     150  case key_gridtype of 
     151    'c':BEGIN 
     152      for i = 0, ny-1, ystride do begin 
     153        plots,  glam[*, i], gphi[*, i], _extra = extra 
     154      endfor 
     155      for i = 0, nx-1, xstride do begin 
     156        plots,  glam[i, *], gphi[i, *], _extra = extra 
     157      endfor 
     158    END 
     159    'e':BEGIN 
     160      shifted = glam[0, 0] LT glam[0, 1] 
     161      glam2 = glam+(glam[1]-glam[0])/2. 
     162      if shifted then begin 
     163        for i = 0, ny-2 do BEGIN 
     164          xx = (transpose([[glam[*, i]], [glam2[*, i]]]))[*] 
     165          yy = (transpose([[gphi[*, i]], [gphi[*, i+1]]]))[*] 
     166          plots, xx[0:2*nx-2], yy[0:2*nx-2], _extra = extra 
     167        ENDFOR 
     168      ENDIF ELSE BEGIN 
     169        for i = 1, ny-1 do BEGIN 
     170          xx = (transpose([[glam[*, i]], [glam2[*, i]]]))[*] 
     171          yy = (transpose([[gphi[*, i]], [gphi[*, i-1]]]))[*] 
     172          plots, xx[0:2*nx-2], yy[0:2*nx-2], _extra = extra 
     173        ENDFOR 
     174      ENDELSE 
     175      for i = 1, (ny-1)/2 do $ 
     176        plots, [glam[0, 2*i-1], glam[0, 2*i]] $ 
     177        , [gphi[0, 2*i-1], gphi[0, 2*i]], _extra = extra 
     178      for i = 0, (ny-2)/2 do $ 
     179        plots, [glam[nx-1, 2*i], glam[nx-1, 2*i+1]] $ 
     180        , [gphi[nx-1, 2*i], gphi[nx-1, 2*i+1]], _extra = extra 
     181    END 
     182  endcase 
    55183 
    56    grille, mask, glam, gphi, gdep, nx, ny 
    57    if (size(mask))[0] EQ 3 then mask = mask[*, *, 0] 
    58    case 1 of 
    59       keyword_set(ocean):BEGIN 
    60          earth = where(mask EQ 0) 
    61          if earth[0] NE -1 then begin 
    62             glam[earth] = !values.f_nan 
    63             gphi[earth] = !values.f_nan 
    64          endif 
    65       END 
    66       keyword_set(earth):BEGIN 
    67          ocean = where(mask EQ 1) 
    68          if ocean[0] NE -1 then begin 
    69             glam[ocean] = !values.f_nan 
    70             gphi[ocean] = !values.f_nan 
    71          endif 
    72       END 
    73       ELSE: 
    74    endcase 
    75    if NOT keyword_set(xstride) then xstride = 1 
    76    if NOT keyword_set(ystride) then ystride = 1 
    77    case key_gridtype of 
    78       'c':BEGIN 
    79          for i = 0, ny-1, ystride do begin 
    80             plots,  glam[*, i], gphi[*, i], color = c_cote, _extra = extra 
    81          endfor 
    82          for i = 0, nx-1, xstride do begin 
    83             plots,  glam[i, *], gphi[i, *], color = c_cote, _extra = extra 
    84          endfor 
    85       END 
    86       'e':BEGIN 
    87          shifted = glam[0, 0] LT glam[0, 1] 
    88          glam2 = glam+(glam[1]-glam[0])/2. 
    89          if shifted then begin 
    90             for i = 0, ny-2 do BEGIN 
    91                xx = (transpose([[glam[*, i]], [glam2[*, i]]]))[*] 
    92                yy = (transpose([[gphi[*, i]], [gphi[*, i+1]]]))[*] 
    93                plots, xx[0:2*nx-2], yy[0:2*nx-2], color = c_cote, _extra = extra 
    94             ENDFOR 
    95          ENDIF ELSE BEGIN 
    96             for i = 1, ny-1 do BEGIN 
    97                xx = (transpose([[glam[*, i]], [glam2[*, i]]]))[*] 
    98                yy = (transpose([[gphi[*, i]], [gphi[*, i-1]]]))[*] 
    99                plots, xx[0:2*nx-2], yy[0:2*nx-2], color = c_cote, _extra = extra 
    100             ENDFOR 
    101          ENDELSE 
    102          for i = 1, (ny-1)/2 do $ 
    103           plots, [glam[0, 2*i-1], glam[0, 2*i]] $ 
    104           , [gphi[0, 2*i-1], gphi[0, 2*i]], color = c_cote, _extra = extra 
    105          for i = 0, (ny-2)/2 do $ 
    106           plots, [glam[nx-1, 2*i], glam[nx-1, 2*i+1]] $ 
    107           , [gphi[nx-1, 2*i], gphi[nx-1, 2*i+1]], color = c_cote, _extra = extra 
    108       END 
    109    endcase 
     184  if keyword_set(key_performance) THEN print, 'temps trace grille', systime(1)-tempsun  
     185  !except = oldexcept 
    110186 
    111    if keyword_set(key_performance) THEN print, 'temps trace grille', systime(1)-tempsun  
    112    return 
     187  return 
    113188end 

  • trunk/ToBeReviewed/GRILLE/u2t.pro

    r12 r13  
    3838;------------------------------------------------------------ 
    3939FUNCTION u2t, temp 
    40 @common 
     40;--------------------------------------------------------- 
     41@cm_4mesh 
     42@cm_4data 
     43@cm_4cal 
     44  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     45@updatenew 
     46  ENDIF 
     47;--------------------------------------------------------- 
    4148   res = temp 
    4249;on force nxt=nxu, etc ... 
    43    premierxt = premierxu 
    44    dernierxt = dernierxu 
    45    premieryt = premieryu 
    46    dernieryt = dernieryu 
     50   firstxt = firstxu 
     51   lastxt = lastxu 
     52   firstyt = firstyu 
     53   lastyt = lastyu 
    4754   nxt = nxu 
    4855   nyt = nyu 
    4956   vargrid = 'T' 
    5057   if NOT keyword_set(valmask) then valmask = 1e20 
    51    lon1 = glamt[premierxt, 0] 
    52    lon2 = glamu[dernierxu, 0] 
     58   lon1 = glamt[firstxt, 0] 
     59   lon2 = glamu[lastxu, 0] 
    5360; 
    5461; cas sur la taille du tableau et application  
    … …  
    6067            taille[1] eq nxu and taille[2] eq nyu: 
    6168            taille[1] eq jpi and taille[2] eq jpj: $ 
    62              res=res[premierxu:dernierxu, premieryu:dernieryu] 
     69             res=res[firstxu:lastxu, firstyu:lastyu] 
    6370            else: $ 
    6471             return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 
    6572         endcase 
    66          mask = (umask())[premierxu:dernierxu, premieryu:dernieryu, 0] 
     73         mask = (umask())[firstxu:lastxu, firstyu:lastyu, 0] 
    6774         terre = where(mask EQ 0)  
    6875         IF terre[0] NE -1 THEN res[terre] = !values.f_nan 
    6976         res = 0.5*(res + shift(res, 1, 0)) 
    70          if NOT (keyword_set(key_periodique) AND nxu EQ jpi) then res[0, *] = !values.f_nan 
    71          mask = tmask[premierxu:dernierxu, premieryu:dernieryu, 0] 
     77         if NOT (keyword_set(key_periodic) AND nxu EQ jpi) then res[0, *] = !values.f_nan 
     78         mask = tmask[firstxu:lastxu, firstyu:lastyu, 0] 
    7279         terre = where(mask EQ 0)  
    7380         IF terre[0] NE -1 THEN res[terre] = valmask 
    … …  
    7784            taille[1] eq nxu and taille[2] eq nyu AND taille[3] EQ nzt: 
    7885            taille[1] eq nxu and taille[2] eq nyu AND taille[3] EQ jpk: $ 
    79              res=res[*, *, premierzt:dernierzt] 
     86             res=res[*, *, firstzt:lastzt] 
    8087            taille[1] eq nxu and taille[2] eq nyu AND taille[3] EQ jpt: 
    8188            taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk: $ 
    82              res=res[premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt] 
     89             res=res[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt] 
    8390            taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpt: $ 
    84              res=res[premierxu:dernierxu, premieryu:dernieryu, *] 
     91             res=res[firstxu:lastxu, firstyu:lastyu, *] 
    8592            else: $ 
    8693             return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 
    8794         ENDCASE 
    8895         if taille[3] EQ jpt then begin 
    89             mask = (umask())[premierxu:dernierxu, premieryu:dernieryu, dernierzt*(nzt NE jpk)] 
     96            mask = (umask())[firstxu:lastxu, firstyu:lastyu, lastzt*(nzt NE jpk)] 
    9097            mask = temporary(mask[*])#replicate(1, jpt) 
    9198            mask = reform(mask, nxu, nyu, jpt, /over) 
    92          ENDIF ELSE mask = (umask())[premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt] 
     99         ENDIF ELSE mask = (umask())[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt] 
    93100         terre = where(temporary(mask) EQ 0)  
    94101         IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan 
    95102         res = 0.5*(res + shift(res, 1, 0, 0)) 
    96          if NOT (keyword_set(key_periodique) AND nxu EQ jpi) then res[0, *, *] = !values.f_nan 
     103         if NOT (keyword_set(key_periodic) AND nxu EQ jpi) then res[0, *, *] = !values.f_nan 
    97104         if taille[3] EQ jpt then BEGIN 
    98             mask = tmask[premierxu:dernierxu, premieryu:dernieryu, dernierzt*(nzt NE jpk)] 
     105            mask = tmask[firstxu:lastxu, firstyu:lastyu, lastzt*(nzt NE jpk)] 
    99106            mask = temporary(mask[*])#replicate(1, jpt) 
    100107            mask = reform(mask, nxu, nyu, jpt, /over) 
    101          ENDIF ELSE mask = tmask[premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt] 
     108         ENDIF ELSE mask = tmask[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt] 
    102109         terre = where(temporary(mask) EQ 0)  
    103110         IF terre[0] NE -1 THEN res[temporary(terre)] = valmask 
    … …  
    107114            taille[1] eq nxu and taille[2] eq nyu AND taille[3] EQ nzt AND taille[4] EQ jpt: 
    108115            taille[1] eq nxu and taille[2] eq nyu AND taille[3] EQ jpk AND taille[4] EQ jpt: $ 
    109              res=res[*, *, premierzt:dernierzt, *] 
     116             res=res[*, *, firstzt:lastzt, *] 
    110117            taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk AND taille[4] EQ jpt: $ 
    111              res=res[premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt, *] 
     118             res=res[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt, *] 
    112119            else: $ 
    113120             return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 
    114121         ENDCASE 
    115          mask = (umask())[premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt] 
     122         mask = (umask())[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt] 
    116123         mask = temporary(mask[*])#replicate(1, jpt) 
    117124         mask = reform(mask, nxu, nyu, nzt, jpt, /over) 
    … …  
    119126         IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan 
    120127         res = 0.5*(res + shift(res, 1, 0, 0, 0)) 
    121          if NOT (keyword_set(key_periodique) AND nxu EQ jpi) then res[0, *, *, *] = !values.f_nan 
    122          mask = tmask[premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt] 
     128         if NOT (keyword_set(key_periodic) AND nxu EQ jpi) then res[0, *, *, *] = !values.f_nan 
     129         mask = tmask[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt] 
    123130         mask = temporary(mask[*])#replicate(1, jpt) 
    124131         mask = reform(mask, nxu, nyu, nzt, jpt, /over) 
    … …  
    128135   endcase 
    129136 
     137  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     138   @updateold 
     139  ENDIF  
     140   
    130141   return, res 
    131142END 

  • trunk/ToBeReviewed/GRILLE/umask.pro

    r12 r13  
    3939;------------------------------------------------------------ 
    4040FUNCTION umask 
    41    tempsun = systime(1)         ; pour key_performance 
    42 @common 
    43    if jpk EQ 1 then begin 
    44       res = tmask*shift(tmask, -1, 0) 
    45       if NOT keyword_set(key_periodique) then res[jpi-1, *] = umaskred 
    46    ENDIF ELSE BEGIN  
    47       res = tmask*shift(tmask, -1, 0, 0) 
    48       if NOT keyword_set(key_periodique) then res[jpi-1, *, *] = umaskred 
    49    ENDELSE 
     41;--------------------------------------------------------- 
     42@cm_4mesh 
     43  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     44@updatenew 
     45  ENDIF 
     46;--------------------------------------------------------- 
     47  tempsun = systime(1)          ; pour key_performance 
    5048; 
    51    if keyword_set(key_performance) THEN print, 'temps umask', systime(1)-tempsun  
     49  CASE size(tmask, /n_dimensions) OF 
     50    2:res = tmask*shift(tmask, -1, 0) 
     51    3:res = tmask*shift(tmask, -1, 0, 0) 
     52  ENDCASE 
     53;    
     54  if NOT keyword_set(key_periodic) then res[jpi-1, *, *] = umaskred 
     55  if keyword_set(key_performance) THEN print, 'temps umask', systime(1)-tempsun  
    5256; 
    53    return, res 
     57  return, res 
    5458end 

  • trunk/ToBeReviewed/GRILLE/v2t.pro

    r12 r13  
    3838;------------------------------------------------------------ 
    3939FUNCTION v2t, temp 
    40 @common 
     40;--------------------------------------------------------- 
     41@cm_4mesh 
     42@cm_4data 
     43@cm_4cal 
     44  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     45@updatenew 
     46  ENDIF 
     47;--------------------------------------------------------- 
    4148   res = temp 
    4249;on force nxt=nxv, etc ... 
    43    premierxt = premierxv 
    44    dernierxt = dernierxv 
    45    premieryt = premieryv 
    46    dernieryt = dernieryv 
     50   firstxt = firstxv 
     51   lastxt = lastxv 
     52   firstyt = firstyv 
     53   lastyt = lastyv 
    4754   nxt = nxv 
    4855   nyt = nyv 
    4956   vargrid = 'T' 
    5057   if NOT keyword_set(valmask) then valmask = 1e20 
    51    lat1 = gphit[0, premieryt] 
    52    lat2 = gphiv[0, dernieryv] 
     58   lat1 = gphit[0, firstyt] 
     59   lat2 = gphiv[0, lastyv] 
    5360 
    5461; cas sur la taille du tableau et application  
    … …  
    6067            taille[1] eq nxv and taille[2] eq nyv: 
    6168            taille[1] eq jpi and taille[2] eq jpj: $ 
    62              res=res[premierxv:dernierxv, premieryv:dernieryv] 
     69             res=res[firstxv:lastxv, firstyv:lastyv] 
    6370            else: $ 
    6471             return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 
    6572         endcase 
    66          mask = (vmask())[premierxv:dernierxv, premieryv:dernieryv, 0] 
     73         mask = (vmask())[firstxv:lastxv, firstyv:lastyv, 0] 
    6774         terre = where(mask EQ 0)  
    6875         IF terre[0] NE -1 THEN res[terre] = !values.f_nan 
    6976         res = 0.5*(res + shift(res, 0, +1)) 
    7077         res[*, 0] = !values.f_nan 
    71          mask = tmask[premierxv:dernierxv, premieryv:dernieryv, 0] 
     78         mask = tmask[firstxv:lastxv, firstyv:lastyv, 0] 
    7279         terre = where(mask EQ 0)  
    7380         IF terre[0] NE -1 THEN res[terre] = valmask 
    … …  
    7784            taille[1] eq nxv and taille[2] eq nyv AND taille[3] EQ nzt: 
    7885            taille[1] eq nxv and taille[2] eq nyv AND taille[3] EQ jpk: $ 
    79              res=res[*, *, premierzt:dernierzt] 
     86             res=res[*, *, firstzt:lastzt] 
    8087            taille[1] eq nxv and taille[2] eq nyv AND taille[3] EQ jpt: 
    8188            taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk: $ 
    82              res=res[premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt] 
     89             res=res[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt] 
    8390            taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpt: $ 
    84              res=res[premierxv:dernierxv, premieryv:dernieryv, *] 
     91             res=res[firstxv:lastxv, firstyv:lastyv, *] 
    8592            else: $ 
    8693             return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 
    8794         ENDCASE 
    8895         if taille[3] EQ jpt then begin 
    89             mask = (vmask())[premierxv:dernierxv, premieryv:dernieryv, dernierzt*(nzt NE jpk)] 
     96            mask = (vmask())[firstxv:lastxv, firstyv:lastyv, lastzt*(nzt NE jpk)] 
    9097            mask = temporary(mask[*])#replicate(1, jpt) 
    9198            mask = reform(mask, nxv, nyv, jpt, /over) 
    92          ENDIF ELSE mask = (vmask())[premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt] 
     99         ENDIF ELSE mask = (vmask())[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt] 
    93100         terre = where(temporary(mask) EQ 0)  
    94101         IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan 
    … …  
    96103         res[*, 0, *] = !values.f_nan 
    97104         if taille[3] EQ jpt then BEGIN 
    98             mask = tmask[premierxv:dernierxv, premieryv:dernieryv, dernierzt*(nzt NE jpk)] 
     105            mask = tmask[firstxv:lastxv, firstyv:lastyv, lastzt*(nzt NE jpk)] 
    99106            mask = temporary(mask[*])#replicate(1, jpt) 
    100107            mask = reform(mask, nxv, nyv, jpt, /over) 
    101          ENDIF ELSE mask = tmask[premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt] 
     108         ENDIF ELSE mask = tmask[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt] 
    102109         terre = where(temporary(mask) EQ 0)  
    103110         IF terre[0] NE -1 THEN res[temporary(terre)] = valmask 
    … …  
    107114            taille[1] eq nxv and taille[2] eq nyv AND taille[3] EQ nzt AND taille[4] EQ jpt: 
    108115            taille[1] eq nxv and taille[2] eq nyv AND taille[3] EQ jpk AND taille[4] EQ jpt: $ 
    109              res=res[*, *, premierzt:dernierzt, *] 
     116             res=res[*, *, firstzt:lastzt, *] 
    110117            taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk AND taille[4] EQ jpt: $ 
    111              res=res[premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt, *] 
     118             res=res[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt, *] 
    112119            else: $ 
    113120             return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 
    114121         ENDCASE 
    115          mask = (vmask())[premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt] 
     122         mask = (vmask())[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt] 
    116123         mask = temporary(mask[*])#replicate(1, jpt) 
    117124         mask = reform(mask, nxv, nyv, nzt, jpt, /over) 
    … …  
    120127         res = 0.5*(res + shift(res, 0, +1, 0, 0)) 
    121128         res[*, 0, *, *] = !values.f_nan 
    122          mask = tmask[premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt] 
     129         mask = tmask[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt] 
    123130         mask = temporary(mask[*])#replicate(1, jpt) 
    124131         mask = reform(mask, nxv, nyv, nzt, jpt, /over) 
    … …  
    128135   endcase 
    129136 
    130    return, res 
     137   IF NOT keyword_set(key_forgetold) THEN BEGIN 
     138   @updateold 
     139   ENDIF  
     140 
     141    return, res 
    131142END 
    132143 

  • trunk/ToBeReviewed/GRILLE/vmask.pro

    r12 r13  
    2727FUNCTION vmask 
    2828@common 
    29    tempsun = systime(1)         ; pour key_performance 
    30    if jpk EQ 1 then begin 
    31       res = tmask*shift(tmask, 0, -1) 
    32       res[*, jpj-1] = vmaskred 
    33    ENDIF ELSE BEGIN  
    34       res = tmask*shift(tmask, 0, -1, 0) 
    35       res[*, jpj-1, *] = vmaskred 
    36    ENDELSE  
    37    if keyword_set(key_performance) THEN print, 'temps vmask', systime(1)-tempsun  
     29  tempsun = systime(1)          ; pour key_performance 
    3830; 
    39    return, res 
     31  CASE size(tmask, /n_dimensions) OF 
     32    2:res = tmask*shift(tmask, 0, -1) 
     33    3:res = tmask*shift(tmask, 0, -1, 0) 
     34  ENDCASE 
     35; 
     36  res[*, jpj-1, *] = vmaskred 
     37  if keyword_set(key_performance) THEN print, 'temps vmask', systime(1)-tempsun  
     38; 
     39  return, res 
    4040end 

  • trunk/ToBeReviewed/GRILLE/whichgrid.pro

    r12 r13  
    44   case filetype of 
    55      'OPA C-Grid.Binary IEEE: Meshmask':BEGIN 
    6          meshlec, filename, /pasblabla, _extra = ex 
     6         meshread, filename, /pasblabla, _extra = ex 
    77      END 
    88      'OPA C-Grid.Net Cdf: Meshmask':BEGIN 
    9          ncdf_meshlec, filename, _extra = ex 
     9         ncdf_meshread, filename, _extra = ex 
    1010      END 
    1111      'Regular 2D.Binary IEEE: mask': 
    … …  
    2828;********************************************************************* 
    2929FUNCTION getgridparameter, top 
    30 @common 
    31  
     30;--------------------------------------------------------- 
     31@cm_4mesh 
     32  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     33@updatenew 
     34  ENDIF 
     35;--------------------------------------------------------- 
    3236   widget_control, widget_info(top, find_by_uname = 'xmesh'), get_value = answer 
    3337   jpiglo = long(answer[0]) 
    … …  
    4852    , get_value = answer 
    4953   key_shift = long(answer[0]) 
    50    widget_control, widget_info(top, find_by_uname = 'key_periodique') $ 
     54   widget_control, widget_info(top, find_by_uname = 'key_periodic') $ 
    5155    , get_value = answer 
    52    key_periodique = long(answer[0]) 
     56   key_periodic = long(answer[0]) 
    5357   widget_control, widget_info(top, find_by_uname = 'triangulation') $ 
    5458    , get_value = answer 
    … …  
    6165          , iyminmesh:iyminmesh, iymaxmesh:iymaxmesh $ 
    6266          , izminmesh:izminmesh, izmaxmesh:izmaxmesh $ 
    63           , key_shift:key_shift, key_periodique:key_periodique $ 
     67          , key_shift:key_shift, key_periodic:key_periodic $ 
    6468          , triangulation:triangulation, boundary:boundary} 
     69 
     70   @updateold 
    6571 
    6672   return, res 
    … …  
    6874;********************************************************************* 
    6975pro showgridparameter, basetop, EDITABLE = editable, _EXTRA = ex 
    70 ; 
    71 @common 
    72 ;------------------------------------------------------------ 
    73 ; 
     76;--------------------------------------------------------- 
     77@cm_4mesh 
     78  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     79@updatenew 
     80  ENDIF 
     81;--------------------------------------------------------- 
    7482   widget_control, basetop, update = 0 
    7583   base=widget_base(basetop, /COLUMN, /align_center, _EXTRA = ex) 
    … …  
    8290   if NOT keyword_set(key_shift) then key_shift = 0 
    8391   nothing = widget_text(basea, value = strtrim(key_shift*(1-keyword_set(editable)),1), uname = 'key_shift', xsize = 4, EDITABLE = editable) 
    84    nothing = widget_label(basea, value = 'key_periodique') 
    85    if NOT keyword_set(key_periodique) then key_periodique = 0 
    86    nothing = widget_text(basea, value = strtrim(key_periodique*(1-keyword_set(editable)),1), uname = 'key_periodique', xsize = 4, EDITABLE = editable) 
     92   nothing = widget_label(basea, value = 'key_periodic') 
     93   if NOT keyword_set(key_periodic) then key_periodic = 0 
     94   nothing = widget_text(basea, value = strtrim(key_periodic*(1-keyword_set(editable)),1), uname = 'key_periodic', xsize = 4, EDITABLE = editable) 
    8795   baseb=widget_base(base, /row, /align_center) 
    8896   nothing = widget_label(baseb, value = 'use a triangulation (y/n) ?') 
    … …  
    159167;********************************************************************* 
    160168PRO whichgrid_event, event 
    161 @common 
     169;--------------------------------------------------------- 
     170@cm_4mesh 
     171  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     172@updatenew 
     173  ENDIF 
     174;--------------------------------------------------------- 
    162175   widget_control, event.id,  get_uvalue = eventuvalue 
    163176   IF chkstru(eventuvalue,'name') EQ 0 THEN return 
    … …  
    181194         if rstrpos(filename, '.pro') EQ strlen(filename)-4 then begin 
    182195            createpro, '@'+strmid(filename, 0, strlen(filename)-4) $ 
    183           , filename = isadirectory(io = homedir, title = 'Bad definition of Homedir') $ 
    184           +'for_createpro.pro' 
     196          , filename = myuniquetmpdir +'for_createpro.pro' 
    185197            showgridparameter, event.handler, group_leader = event.handler,/frame, uname = 'showgridparameter' 
    186198         ENDIF ELSE BEGIN 
    … …  
    204216         if rstrpos(filename, '.pro') EQ strlen(filename)-4 then begin 
    205217            createpro, '@'+strmid(filename, 0, strlen(filename)-4) $ 
    206           , filename = isadirectory(io = homedir, title = 'Bad definition of Homedir') $ 
    207           +'for_createpro.pro' 
     218          , filename = myuniquetmpdir +'for_createpro.pro' 
    208219            showgridparameter, event.handler, group_leader = event.handler,/frame, uname = 'showgridparameter' 
    209220         ENDIF ELSE BEGIN 
    … …  
    223234         case filetype of 
    224235            'OPA C-Grid.Binary IEEE: Meshmask':BEGIN 
    225                meshlec, filename, /pasblabla, /getdimensions 
     236               meshread, filename, /pasblabla, /getdimensions 
    226237               showgridparameter, event.handler, /editable, group_leader = event.handler,/frame, uname = 'showgridparameter' 
    227238            END 
    228239            'OPA C-Grid.Net Cdf: Meshmask':BEGIN 
    229                ncdf_meshlec, filename, /getdimensions 
     240               ncdf_meshread, filename, /getdimensions 
    230241               showgridparameter, event.handler, /editable, group_leader = event.handler,/frame, uname = 'showgridparameter' 
    231242            END 
    … …  
    269280END 
    270281;********************************************************************* 
    271 FUNCTION whichgrid, name, IODIRECTORY = iodirectory, PARENT = parent, _EXTRA = ex 
    272 @common 
    273 ; 
     282FUNCTION whichgrid, name, PARENT = parent, _EXTRA = ex 
     283;--------------------------------------------------------- 
     284@cm_4mesh 
     285  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     286@updatenew 
     287  ENDIF 
     288;--------------------------------------------------------- 
    274289; 
    275290   if n_elements(name) NE 0 then begin 
    276       filename = isafile(filename = name) 
     291      filename = isafile(filename = name, IODIRECTORY = iodir, _extra = ex) 
    277292      if size(filename, /type) NE 7 then return, -1 
    278293   ENDIF ELSE filaname = 'no file' 
    … …  
    301316      if rstrpos(filename, '.pro') EQ strlen(filename)-4 then begin 
    302317         createpro, '@'+strmid(filename, 0, strlen(filename)-4) $ 
    303           , filename = isadirectory(io = homedir, title = 'Bad definition of Homedir') $ 
    304           +'for_createpro.pro' 
     318          , filename = myuniquetmpdir +'for_createpro.pro' 
    305319         showgridparameter, base, group_leader = base,/frame, uname = 'showgridparameter' 
    306320      ENDIF ELSE BEGIN 
Note: See TracChangeset for help on using the changeset viewer.