Ignore:
Timestamp:
05/02/06 15:32:01 (18 years ago)
Author:
pinsard
Message:

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

File:
1 copied

Legend:

Unmodified
Added
Removed

  • trunk/ToBeReviewed/TRIANGULATION/section.pro

    r27 r29  
    3232;------------------------------------------------------------ 
    3333 
    34 PRO section, field, res, glamaxe, gphiaxe, ENDPOINTS = endpoints, BOITE = boite, TYPE = type, DIREC = direc 
    35 @common 
    36 ;  
    37 ;--------------------- 
    38    array = litchamp(field) 
    39 ;--------------------- 
    40 ; definition de boite en fonction de endpoints 
     34PRO section, field, res, glamaxe, gphiaxe, ENDPOINTS = endpoints $ 
     35             , BOXZOOM = boxzoom, TYPE = type, WDEPTH = wdepth $ 
     36             , DIREC = direc, SHOWBUILD = showbuild, ONLYBOX = onlybox $ 
     37             , _extra = ex 
     38; 
     39;--------------------------------------------------------- 
     40; include common 
     41@cm_4mesh 
     42@cm_4data 
     43@cm_4cal 
     44  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     45@updatenew 
     46@updatekwd 
     47  ENDIF 
     48;-------------------------------------------------------------- 
     49;--------------------- 
     50;--------------------- 
     51; definition de boxzoom en fonction de endpoints 
    4152; puis redefinition du domaine 
    42    boite2d = [min([endpoints[0], endpoints[2]]), max([endpoints[0], endpoints[2]]), min([endpoints[1], endpoints[3]]), max([endpoints[1], endpoints[3]])] 
    43 ; 
    44    minprof = 0 
    45    profdefault = 200 
    46    if n_elements(type) EQ 0 then type = 'nothing' 
    47    Case N_Elements(Boite) OF 
    48       1:boite=[boite2d, minprof, boite[0]] 
    49       2:boite=[boite2d, boite[0],boite[1]] 
    50       Else:$  
    51        if strpos(type, 'z') NE -1 THEN boite=[boite2d, minprof, profdefault] $ 
    52       ELSE boite=[boite2d, prof1, prof2] 
    53    ENDCASE 
    54    if strpos(type, 'z') NE -1 THEN BEGIN 
    55       !y.range = [boite[5], boite[4]] ;on garde les yranges (axe z) avant de changer la boite. 
    56       profmax = boite[5] 
    57       vraiboite = boite 
    58       if vargrid EQ 'W' then gdep = gdepw ELSE gdep = gdept 
    59 ; check with vertical grid limits (nearest level) 
    60       gwork = gdep 
    61 ; check the increse or decrese of the z axis 
    62       IF gwork[1] LE gwork[0] THEN gwork = reverse(gdep, 1) 
    63       niveauprof = where(gwork ge boite[5]) & niveauprof = niveauprof[0] 
    64       if niveauprof NE -1 then boite[5] = gwork[niveauprof]+1 
    65    ENDIF 
    66    domdef, boite, GRILLE=[vargrid], /findalways, _extra = ex 
     53  boxzoom2d = [min([endpoints[0], endpoints[2]], max = ma02), ma02 $ 
     54               , min([endpoints[1], endpoints[3]], max = ma13), ma13] 
     55; 
     56  minprof = 0. 
     57  profdefault = 200. 
     58  if n_elements(type) EQ 0 then type = 'nothing' 
     59  Case N_Elements(Boxzoom) OF 
     60    0:localbox = [boxzoom2d, minprof, profdefault] 
     61    1:localbox = [boxzoom2d, minprof, boxzoom[0]] 
     62    2:localbox = [boxzoom2d, boxzoom[0]] 
     63    4:if strpos(type, 'z') NE -1 THEN $ 
     64      localbox = [boxzoom2d, minprof, profdefault] ELSE localbox = boxzoom2d 
     65    5:localbox = [boxzoom2d, minprof, boxzoom[4]] 
     66    6:localbox = [boxzoom2d, boxzoom[4:5]] 
     67    Else:BEGIN 
     68      print, report('Bad definition of the box') 
     69      stop 
     70    END 
     71  ENDCASE 
     72  nelbox = n_elements(localbox) 
     73; 
     74  if keyword_set(wdepth) then grillechoice = [vargrid, 'W'] $ 
     75  ELSE grillechoice = vargrid 
     76  domdef, localbox, GRIDTYPE = grillechoice, /findalways, _extra = ex 
     77  grille, -1, -1, -1, -1, nx, ny 
     78; if less than 10 points where found, we apply domdef over the whole domain  
     79; -> problem... why 10 points as a test value??? 
     80; how can we find a good test value? 
     81  IF nx * ny LE 10 THEN domdef, GRIDTYPE = grillechoice, _extra = ex 
    6782; on redefinit lon1, ... au cas ou findalways ait ete utilise ds domdef 
    68    lon1 = min([endpoints[0], endpoints[2]]) 
    69    lon2 = max([endpoints[0], endpoints[2]]) 
    70    lat1 = min([endpoints[1], endpoints[3]]) 
    71    lat2 = max([endpoints[1], endpoints[3]]) 
    72 ; on recupere la grille sur la boite 
    73    grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz,dernierx, derniery, dernierz 
     83  lon1 = min([endpoints[0], endpoints[2]], max = lon2) 
     84  lat1 = min([endpoints[1], endpoints[3]], max = lat2) 
     85; we extend the box along the z axis -> i that way the plot will be drawn 
     86; until its bottom part. 
     87  if strpos(type, 'z') NE -1 THEN BEGIN 
     88;on garde les yranges (axe z) avant de changer la boxzoom. 
     89    !y.range = [localbox[nelbox-1], localbox[nelbox-2]]  
     90    if vargrid EQ 'W' OR keyword_set(wdepth) then BEGIN 
     91      firstzw = 0 > (firstzw-1) 
     92      lastzw = (lastzw+1) < (jpk-1) 
     93      nzw = lastzw - firstzw + 1 
     94    ENDIF ELSE BEGIN 
     95      firstzt = 0 > (firstzt-1) 
     96      lastzt = (lastzt+1) < (jpk-1) 
     97      nzt = lastzt - firstzt + 1 
     98    ENDELSE 
     99    IF NOT keyword_set(key_forgetold) THEN BEGIN 
     100     @updateold 
     101   ENDIF  
     102  ENDIF 
     103  !y.range = [localbox[nelbox-1], localbox[nelbox-2]]  
     104; on recupere la grille sur la boxzoom 
     105; 
     106  IF NOT keyword_set(onlybox) THEN saveboxparam, 'boxparam4section.dat' 
     107  grille, -1, -1, -1, -1, nx, ny, nz $ 
     108          , firstx, firsty, firstz, lastx, lasty, lastz 
     109; the extend the box in the x and y direction in order to maximise 
     110; the number of triangles used to build the section 
     111  firstx = 0 > (firstx - 1) 
     112  lastx = (lastx + 1) < (jpi -1) 
     113  firsty = 0 > (firsty - 1) 
     114  lasty = (lasty + 1) < (jpj -1) 
     115   
     116  domdef, firstx, lastx, firsty, lasty, firstz, lastz $ 
     117          , /index, gridtype = vargrid 
     118 
     119  IF keyword_set(onlybox) THEN return 
     120; 
     121  grille, mask, glam, gphi, gdep, nx, ny, nz $ 
     122          , firstx, firsty, firstz, lastx, lasty, lastz 
     123 
    74124;--------------------- 
    75125; on definit la triangulation qui va nous permetre de determiner la 
    … …  
    77127; aussi bien que sur la mer. suivant le sens de la section -plutot 
    78128; longitude ou plutot latitude- on definit la facon de trianguler 
    79    if strpos(type, 'x') NE -1 then BEGIN  
    80       downward = (lindgen(nx, ny))[0:nx-2, 0:ny-2] 
    81       tri = definetri(nx, ny, (downward)[*])  
    82    ENDIF ELSE tri = definetri(nx, ny) 
     129  if strpos(type, 'x') NE -1 then BEGIN  
     130    downward = (lindgen(nx, ny))[0:nx-2, 0:ny-2] 
     131    tri = definetri(nx, ny, (downward)[*])  
     132  ENDIF ELSE tri = definetri(nx, ny) 
     133; If we have an irregular grid that is periodic, then it is possible that 
     134; some of the triangle have a very large size (neighborg points on the 
     135; sphere but far away when doing the projection) and should not be 
     136; taken into account..  
     137  IF keyword_set(key_irregular) AND keyword_set(key_periodic) THEN BEGIN 
     138    glamtri = glam[tri] 
     139    glamtri = abs(glamtri - shift(glamtri, 1, 0)) 
     140    good = temporary(glamtri) LT (10.*max(glam)/nx) 
     141    good = where(total(temporary(good), 1) EQ 3) 
     142    tri = (temporary(tri))[*, temporary(good)] 
     143  ENDIF 
    83144;--------------------- 
    84145; equation de la droite suivant laquelle on fait la section 
    85146;--------------------- 
    86    abc = linearequation(endpoints[0:1], endpoints[2:3]) 
    87    glamtri = glam[tri] 
    88    gphitri = gphi[tri] 
     147  abc = linearequation(endpoints[0:1], endpoints[2:3]) 
     148  glamtri = glam[tri] 
     149  gphitri = gphi[tri] 
    89150; quels sont les points de la triangulation qui sont au dessus et au 
    90151; dessous de la droite ? 
    91    if abc[1] NE 0 THEN $ 
     152  if abc[1] NE 0 THEN $ 
    92153    test = temporary(gphitri) GE (-abc[0]/abc[1]*temporary(glamtri)-abc[2]/abc[1]) $ 
    93    ELSE test = temporary(glamtri) GE (-abc[1]/abc[0]*temporary(gphitri)-abc[2]/abc[0]) 
    94  
    95    zero123 = total(test, 1) 
     154  ELSE test = temporary(glamtri) GE (-abc[1]/abc[0]*temporary(gphitri)-abc[2]/abc[0]) 
     155 
     156  zero123 = total(test, 1) 
    96157; to keep: triangles de la triangulation qui sont a cheval sur la droite. 
    97    tokeep1 = where(zero123 EQ 1) 
    98    tokeep2 = where(temporary(zero123) EQ 2) 
    99    tokeep = [tokeep1, tokeep2] 
    100  
    101    test = test[*, tokeep] 
    102    tri = tri[*, tokeep] 
     158  tokeep1 = where(zero123 EQ 1) 
     159  tokeep2 = where(temporary(zero123) EQ 2) 
     160  tokeep = [tokeep1, tokeep2] 
     161 
     162  test = test[*, tokeep] 
     163  tri = tri[*, tokeep] 
    103164; quel est le sommet du triangle qui est seul d''un cote de la droite? 
    104    single1 = where(test[*, 0:n_elements(tokeep1)-1] EQ 1) 
    105    single1 = single1-(single1/3)*3 
    106    single2 = where(test[*, n_elements(tokeep1):n_elements(tokeep)-1] EQ 0) 
    107    single2 = single2-(single2/3)*3 
    108  
    109    undefine, tokeep 
    110    undefine, tokeep1 
    111    undefine, tokeep2 
    112    undefine, test 
    113  
    114    single = [temporary(single1), temporary(single2)] 
     165  single1 = where(test[*, 0:n_elements(tokeep1)-1] EQ 1) 
     166  single1 = single1-(single1/3)*3 
     167  single2 = where(test[*, n_elements(tokeep1):n_elements(tokeep)-1] EQ 0) 
     168  single2 = single2-(single2/3)*3 
     169 
     170  undefine, tokeep 
     171  undefine, tokeep1 
     172  undefine, tokeep2 
     173  undefine, test 
     174 
     175  single = [temporary(single1), temporary(single2)] 
    115176; points1 le point du triangles qui est seul d''un cote de la droite. 
    116177; point2 l''autre point du triangle de l''autre cote de la droite 
    117    point1 = [single, single] 
    118    point2 = [single EQ 0, 1 + (single LE 1)] 
    119  
    120    undefine,  single 
    121  
    122    ntri = (size(tri))[2] 
    123    index = [lindgen(ntri), lindgen(ntri)] 
    124  
    125    points1 = tri[point1, index] 
    126    points2 = tri[point2, temporary(index)] 
     178  point1 = [single, single] 
     179  point2 = [single EQ 0, 1 + (single LE 1)] 
     180 
     181  undefine,  single 
     182 
     183  ntri = (size(tri))[2] 
     184  index = [lindgen(ntri), lindgen(ntri)] 
     185 
     186  points1 = tri[point1, index] 
     187  points2 = tri[point2, temporary(index)] 
    127188; points : complexe contenant les couples de points de part et 
    128189; d''autre de la droite. Ils faut supprimer les doublons. 
    129    points = dcomplex(points1, points2) 
    130    points = points[uniq(points, sort(points))] 
    131    symetrique = dcomplex(imaginary(points), double(points)) 
    132    points = points[where(points-shift(temporary(symetrique), 1) NE 0)] 
     190  points = dcomplex(points1, points2) 
     191  points = points[uniq(points, sort(points))] 
     192  symetrique = dcomplex(imaginary(points), double(points)) 
     193  points = points[where(points-shift(temporary(symetrique), 1) NE 0)] 
    133194; points1 les coordnnees du point du triangles qui est seul d''un cote de la droite. 
    134195; point2 les coordnnees de l''autre point du triangle de l''autre cote de la droite 
    135    points1 = complex(glam[   double(points)], gphi[   double(points)]) 
    136    points2 = complex(glam[imaginary(points)], gphi[imaginary(points)]) 
     196  points1 = complex(glam[   double(points)], gphi[   double(points)]) 
     197  points2 = complex(glam[imaginary(points)], gphi[imaginary(points)]) 
    137198; droites les equations des droites dont on cherche l''intersection 
    138199; avec la section. 
    139    droites = linearequation(points1, points2) 
    140    inter = lineintersection(droites, abc#replicate(1, n_elements(points1) ) ) 
     200  droites = linearequation(points1, points2) 
     201  inter = lineintersection(droites, abc#replicate(1, n_elements(points1) ) ) 
    141202 
    142203; les ccordonnes geographiques des points que l''on cherche sur la section. 
    143    glamaxe = float(inter) 
    144    gphiaxe = imaginary(inter) 
     204  glamaxe = float(inter) 
     205  gphiaxe = imaginary(inter) 
    145206; on les range ds l''ordre croissant entre les bornes de la section 
    146    if strpos(type, 'x') NE -1 then BEGIN  
    147       sort = sort(glamaxe) 
    148       glamaxe = glamaxe[sort] 
    149       inbox = where(glamaxe GE lon1 AND glamaxe LE lon2) 
    150       glamaxe = glamaxe[inbox] 
    151       sort = sort[inbox] 
    152       gphiaxe = gphiaxe[sort] 
    153    ENDIF ELSE BEGIN 
    154       sort = sort(gphiaxe) 
    155       gphiaxe = gphiaxe[sort] 
    156       inbox = where(gphiaxe GE lat1 AND gphiaxe LE lat2) 
    157       gphiaxe = gphiaxe[inbox] 
    158       sort = sort[inbox] 
    159       glamaxe = glamaxe[sort] 
    160    ENDELSE 
    161    points = points[sort] 
    162    points1 = points1[sort] 
    163    points2 = points2[sort] 
    164    inter = inter[sort] 
    165    poids = abs(points2-inter)/abs(points2-points1) 
    166  
    167 ;--------------------- 
    168    array = fitintobox(array) 
    169    if array[0] EQ -1 THEN BEGIN 
    170       res = -1 
    171       return 
    172    ENDIF 
    173    if n_elements(valmask) EQ 0 THEN valmask = 1e20 
    174    taille=size(array) 
    175    if jpt GT 1 AND taille[0] GE 3 AND strpos(type, 't') EQ -1 then BEGIN 
    176       direc = 't' 
    177       array = grossemoyenne(array, 't') 
    178       taille=size(array) 
    179       jpt = 1 
    180    ENDIF 
    181    case 1 of 
     207  if strpos(type, 'x') NE -1 then BEGIN  
     208    sort = sort(glamaxe) 
     209    glamaxe = glamaxe[sort] 
     210    inbox = where(glamaxe GE lon1 AND glamaxe LE lon2) 
     211    glamaxe = glamaxe[inbox] 
     212    sort = sort[inbox] 
     213    gphiaxe = gphiaxe[sort] 
     214  ENDIF ELSE BEGIN 
     215    sort = sort(gphiaxe) 
     216    gphiaxe = gphiaxe[sort] 
     217    inbox = where(gphiaxe GE lat1 AND gphiaxe LE lat2) 
     218    gphiaxe = gphiaxe[inbox] 
     219    sort = sort[inbox] 
     220    glamaxe = glamaxe[sort] 
     221  ENDELSE 
     222  points = points[sort] 
     223  points1 = points1[sort] 
     224  points2 = points2[sort] 
     225  inter = inter[sort] 
     226  poids = abs(points2-inter)/abs(points2-points1) 
     227 
     228;--------------------- 
     229  array = litchamp(field) 
     230  array = fitintobox(array) 
     231  if array[0] EQ -1 THEN BEGIN 
     232    res = -1 
     233    return 
     234  ENDIF 
     235  if n_elements(valmask) EQ 0 THEN valmask = 1e20 
     236  taille = size(array) 
     237  if jpt GT 1 AND taille[0] GE 3 AND strpos(type, 't') EQ -1 then BEGIN 
     238    direc = 't' 
     239    array = grossemoyenne(array, 't') 
     240    taille = size(array) 
     241    jpt = 1 
     242  ENDIF 
     243  case 1 of 
    182244;---------------------------------------------------------------------------- 
    183245;xy 
    184246;---------------------------------------------------------------------------- 
    185       taille[0] EQ 2:BEGIN 
    186          value1 = array[double(points)] 
    187          terre = where(value1 GT valmask/10) 
    188          if terre[0] NE -1 then value1[terre] = !values.f_nan 
    189          value2 = array[imaginary(points)] 
    190          terre = where(value2 GT valmask/10) 
    191          if terre[0] NE -1 then value2[terre] = !values.f_nan 
    192          res = poids*value1+(1-poids)*value2 
    193       END 
     247    taille[0] EQ 2:BEGIN 
     248      value1 = array[double(points)] 
     249      terre = where(value1 GT valmask/10) 
     250      if terre[0] NE -1 then value1[terre] = !values.f_nan 
     251      value2 = array[imaginary(points)] 
     252      terre = where(value2 GT valmask/10) 
     253      if terre[0] NE -1 then value2[terre] = !values.f_nan 
     254      res = poids*value1+(1-poids)*value2 
     255    END 
    194256;---------------------------------------------------------------------------- 
    195257;xyz 
    196258;---------------------------------------------------------------------------- 
    197       taille[0] EQ 3 AND jpt EQ 1:BEGIN 
    198          npoints = n_elements(points)  
    199          index = double(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz) 
    200          value1 = array[index] 
    201          terre = where(value1 GT valmask/10) 
    202          if terre[0] NE -1 then value1[terre] = !values.f_nan 
    203          index = imaginary(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz) 
    204          value2 = array[index] 
    205          terre = where(value2 GT valmask/10) 
    206          if terre[0] NE -1 then value2[terre] = !values.f_nan 
    207          poids = poids#replicate(1, nz) 
    208          res = poids*value1+(1-poids)*value2 
     259    taille[0] EQ 3 AND jpt EQ 1:BEGIN 
     260      npoints = n_elements(points)  
     261      index = double(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz) 
     262      value1 = array[index] 
     263      terre = where(value1 GT valmask/10) 
     264      if terre[0] NE -1 then value1[terre] = !values.f_nan 
     265      index = imaginary(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz) 
     266      value2 = array[index] 
     267      terre = where(value2 GT valmask/10) 
     268      if terre[0] NE -1 then value2[terre] = !values.f_nan 
     269      poids = poids#replicate(1, nz) 
     270      res = poids*value1+(1-poids)*value2 
    209271; moyenne suivant z ? 
    210          if strpos(type, 'z') EQ -1 then begin 
    211             nan = where(finite(res) EQ 0) 
    212             if vargrid EQ 'W' then e3 = e3w[premierzw:dernierzw] ELSE e3 = e3t[premierzt:dernierzt] 
    213             weight = replicate(1, npoints)#e3 
    214             if nan[0] NE -1 then weight[nan] = !values.f_nan 
    215             totalweight = total(weight, 2, /nan) 
    216             zero = where(totalweight EQ 0) 
    217             if zero[0] NE -1 then totalweight[zero] = !values.f_nan 
    218             res = total(res*weight, 2, /nan)/totalweight 
    219             direc = 'z'+string(byte(testvar(var=toto))) 
    220          endif 
    221       END 
     272      if strpos(type, 'z') EQ -1 then begin 
     273        nan = where(finite(res) EQ 0) 
     274        if vargrid EQ 'W' then e3 = e3w[firstzw:lastzw] ELSE e3 = e3t[firstzt:lastzt] 
     275        weight = replicate(1, npoints)#e3 
     276        if nan[0] NE -1 then weight[nan] = !values.f_nan 
     277        totalweight = total(weight, 2, /nan) 
     278        zero = where(totalweight EQ 0) 
     279        if zero[0] NE -1 then totalweight[zero] = !values.f_nan 
     280        res = total(res*weight, 2, /nan)/totalweight 
     281        direc = 'z'+string(byte(testvar(var = toto))) 
     282      endif 
     283    END 
    222284;---------------------------------------------------------------------------- 
    223285;xyt 
    224286;---------------------------------------------------------------------------- 
    225       taille[0] EQ 3 AND jpt NE 1:BEGIN 
    226          npoints = n_elements(points)  
    227          index = double(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt) 
    228          value1 = array[index] 
    229          terre = where(value1 GT valmask/10) 
    230          if terre[0] NE -1 then value1[terre] = !values.f_nan 
    231          index = imaginary(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt) 
    232          value2 = array[index] 
    233          terre = where(value2 GT valmask/10) 
    234          if terre[0] NE -1 then value2[terre] = !values.f_nan 
    235          poids = poids#replicate(1, jpt) 
    236          res = poids*value1+(1-poids)*value2 
    237       END 
     287    taille[0] EQ 3 AND jpt NE 1:BEGIN 
     288      npoints = n_elements(points)  
     289      index = double(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt) 
     290      value1 = array[index] 
     291      terre = where(value1 GT valmask/10) 
     292      if terre[0] NE -1 then value1[terre] = !values.f_nan 
     293      index = imaginary(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt) 
     294      value2 = array[index] 
     295      terre = where(value2 GT valmask/10) 
     296      if terre[0] NE -1 then value2[terre] = !values.f_nan 
     297      poids = poids#replicate(1, jpt) 
     298      res = poids*value1+(1-poids)*value2 
     299    END 
    238300;---------------------------------------------------------------------------- 
    239301;xyzt 
    240302;---------------------------------------------------------------------------- 
    241       taille[0] EQ 4:BEGIN 
    242          npoints = n_elements(points)  
    243          index = double(points)#replicate(1, nz*jpt)+replicate(nx*ny, npoints)#lindgen(nz*jpt) 
    244          index = reform(index, npoints, nz, jpt, /over) 
    245          value1 = array[index] 
    246          terre = where(value1 GT valmask/10) 
    247          if terre[0] NE -1 then value1[terre] = !values.f_nan 
    248          index = imaginary(points)#replicate(1, nz*jpt)+replicate(nx*ny, npoints)#lindgen(nz*jpt) 
    249          index = reform(index, npoints, nz, jpt, /over) 
    250          value2 = array[index] 
    251          terre = where(value2 GT valmask/10) 
    252          if terre[0] NE -1 then value2[terre] = !values.f_nan 
    253          poids = poids#replicate(1, nz*jpt) 
    254          poids = reform(poids, npoints, nz, jpt, /over) 
    255          res = poids*value1+(1-poids)*value2 
     303    taille[0] EQ 4:BEGIN 
     304      npoints = n_elements(points)  
     305      index = double(points)#replicate(1, nz*jpt)+replicate(nx*ny, npoints)#lindgen(nz*jpt) 
     306      index = reform(index, npoints, nz, jpt, /over) 
     307      value1 = array[index] 
     308      terre = where(value1 GT valmask/10) 
     309      if terre[0] NE -1 then value1[terre] = !values.f_nan 
     310      index = imaginary(points)#replicate(1, nz*jpt)+replicate(nx*ny, npoints)#lindgen(nz*jpt) 
     311      index = reform(index, npoints, nz, jpt, /over) 
     312      value2 = array[index] 
     313      terre = where(value2 GT valmask/10) 
     314      if terre[0] NE -1 then value2[terre] = !values.f_nan 
     315      poids = poids#replicate(1, nz*jpt) 
     316      poids = reform(poids, npoints, nz, jpt, /over) 
     317      res = poids*value1+(1-poids)*value2 
    256318; moyenne suivant z ? 
    257          if strpos(type, 'z') EQ -1 then begin 
    258             nan = where(finite(res) EQ 0) 
    259             if vargrid EQ 'W' then e3 = e3w[premierzw:dernierzw] ELSE e3 = e3t[premierzt:dernierzt] 
    260             weight = replicate(1, npoints)#e3 
    261             weight = weight[*]#replicate(1, jpt) 
    262             weight = reform(weight, npoints, nz, jpt, /over) 
    263             if nan[0] NE -1 then weight[nan] = !values.f_nan 
    264             totalweight = total(weight, 2, /nan) 
    265             zero = where(totalweight EQ 0) 
    266             if zero[0] NE -1 then totalweight[zero] = !values.f_nan 
    267             res = total(res*weight, 2, /nan)/totalweight 
    268             direc = 'z'+string(byte(testvar(var=toto))) 
    269          endif 
    270       END 
    271    endcase 
    272 ;--------------------- 
    273  
    274    terre = where(finite(res) EQ 0) 
    275    if terre[0] NE -1 then res[terre] = valmask 
    276  
    277  
    278 ;    plt, tmask[*,*,0],/nocouleur, /rempli, title = '', subtitle = '' 
    279 ;    !p.title='' 
    280 ;    !p.subtitle='' 
    281 ;    dessinetri 
    282 ;    plot,[endpoints[0],endpoints[2]],[endpoints[1],endpoints[3]],/noerase,color=50 
    283  
    284 ;    plot,float(points1),imaginary(points1),/noerase,color=150,psym=1 
    285 ;    plot,float(points2),imaginary(points2),/noerase,color=150,psym=1 
    286 ;    plot,float(inter),imaginary(inter),/noerase,color=250,psym=1 
    287  
    288 ;    plot,[float(points1[0])],[imaginary(points1[0])],/noerase,color=150,psym=1 
    289 ;    plot,[float(points2[0])],[imaginary(points2[0])],/noerase,color=150,psym=1 
    290 ;    plot,[float(inter[0])],[imaginary(inter[0])],/noerase,color=250,psym=1 
    291  
    292  
    293 ;--------------------- 
    294    return 
     319      if strpos(type, 'z') EQ -1 then begin 
     320        nan = where(finite(res) EQ 0) 
     321        if vargrid EQ 'W' then e3 = e3w[firstzw:lastzw] ELSE e3 = e3t[firstzt:lastzt] 
     322        weight = replicate(1, npoints)#e3 
     323        weight = weight[*]#replicate(1, jpt) 
     324        weight = reform(weight, npoints, nz, jpt, /over) 
     325        if nan[0] NE -1 then weight[nan] = !values.f_nan 
     326        totalweight = total(weight, 2, /nan) 
     327        zero = where(totalweight EQ 0) 
     328        if zero[0] NE -1 then totalweight[zero] = !values.f_nan 
     329        res = total(res*weight, 2, /nan)/totalweight 
     330        direc = 'z'+string(byte(testvar(var = toto))) 
     331      endif 
     332    END 
     333  endcase 
     334;--------------------- 
     335 
     336  terre = where(finite(res) EQ 0) 
     337  if terre[0] NE -1 then res[terre] = valmask 
     338 
     339  if n_elements(showbuild) then BEGIN  
     340    winsave = !window 
     341    psave = !p 
     342    xsave = !x 
     343    ysave = !y 
     344    plt, findgen(nx, ny), /nodata, /nofill, /rempli, title = '', subtitle = '', coast_thick = 2, window = showbuild 
     345    !p.title = '' 
     346    !p.subtitle = '' 
     347 
     348    plots, [endpoints[0], endpoints[2]], [endpoints[1], endpoints[3]], color = 50 
     349    plots, [endpoints[0], endpoints[2]], [endpoints[1], endpoints[3]], color = 50, psym = 2, thick = 2 
     350 
     351    FOR i = 0, n_elements(points1)-1 DO $ 
     352      plots, [float(points1[i]), float(points2[i])] $ 
     353             , [imaginary(points1[i]), imaginary(points2[i])], color = 150 
     354 
     355    plots, float(points1), imaginary(points1), color = 150, psym = 1 
     356    plots, float(points2), imaginary(points2), color = 150, psym = 1 
     357    plots, float(inter), imaginary(inter), color = 250, psym = 1 
     358     
     359    IF terre[0] NE -1 THEN plots, float(inter[terre]), imaginary(inter[terre]), color = 0, psym = 1 
     360     
     361;      dummy = '' 
     362;      read, dummy,  prompt = 'press return to continue' 
     363    IF !d.name EQ 'PS' THEN erase ELSE wset, winsave 
     364    !p = psave 
     365    !x = xsave 
     366    !y = ysave 
     367  ENDIF  
     368 
     369  restoreboxparam, 'boxparam4section.dat' 
     370 
     371;--------------------- 
     372  return 
    295373end 
Note: See TracChangeset for help on using the changeset viewer.